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ABSTRACT 


A  nonlinear  optimal  feedback  control  scheme  for  controlling  a  vehicle 
re-entering  the  earth's  atmosphere  from  lunar  return  initial  conditions 
is  reported. 

The  optimal  feedback  control  law  used  in  the  scheme  is  obtained  from  a  multi 
dimensional  surface  fit  of  the  control  function  for  several  optimal  trajectories. 
Partials  of  the  control  with  respect  to  the  state  vector  are  included  in  the  fit¬ 
ting  procedure.  The  functional  minimized  by  the  trajectories  is  the  total 
(convective  plus  radiative)  stagnation  point  heat. 

The  feedback  control  scheme  is  developed,  and  several  re-entry  trajectories 
are  simulated.  Modest  increases  in  total  heat  from  optimal  values  are  ob¬ 
served,  and  large  (although  tolerable)  terminal  point  errors  occur.  It  is 
believed  that  the  terminal  errors  can  be  greatly  reduced,  if  necessary. 

A  powerful  predictor  scheme  is  developed  which  allows  optimal  trajectories 
to  be  changed  as  a  function  of  a  parameter.  This  is  used  to  extend  the  range 
of  an  optimal  trajectory,  to  perform  an  "absolute  minimum"  test,  and  to  map 
the  optimal  re-entry  corridor. 


Sufficiency  tests  for  a  relative  minimum  are  mechanized,  and  it  is  shown  that 
the  trajectories  considered  are  minimizing  paths.  The  optimization  method 
is  extended  to  inciude  the  bounded  state-coordinate  problem. 
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SECTION  I 
INTRODUCTION 


Many  schemes  have  appeared,  over  the  past  few  years,  for  tne  control  of 
vehicles  re-entering  the  earth's  atmosphere.  They  range  over  a  wide  spec¬ 
trum  of  missions,  from  the  attainment  of  a  given  landing  site  to  rather 
complex  in-flight  maneuvers,  and,  in  operation,  may  be  completely  auto¬ 
matic  or  may  require  a  pilot  in  the  loop.  The  schemes  have  been  developed 
for  a  variety  of  vehicles,  onboard  sensors,  and  onboard  computing  capabilities 
and,  additionally,  cover  a  wide  range  of  initial  re-entry  conditions.  A  good 
survey  is  given  by  Wingrove  in  Reference  15. 

Under  the  heading  of  optimal  control,  re-entry  schemes  become. sparse.  The 
only  seemingly  applicable  approach  published  to  date  is  the  "neighboring  op¬ 
timum  control  scheme"  of  Reference  8.  This  is  a  linear  control  scheme  based 
on  small  perturbations  from  a  nominal  optimal  trajectory.  Previous  studies 
(Reference  1)  indicate  that  such  a  method  may  not  be  applicable  for  re-entry 
control,  although  for  other  applications  it  may  be  perfectly  suitable.  The 
objective  of  this  report  is  to  describe  a  nonlinear  optimal  feedback  control 
scheme,  developed  and  simulated  during  the  contract  period. 

Figure  1-i  illustrates  the  approach.  The  vehicle  is  the  plant  to  be  controlled 
through  application  of  the  control  vector  u(t),  ana  sensors  measure  parameters 
of  the  motion  a.  The  navigator  box  supplies  the  siaie  vecior  and  its  time 
derivative  at  discrete  and  equally  spaced  instants  of  time.  (The  time  incre¬ 
ment  may  be  made  variable  if  necessary. )  A  predictor  equation  estimates 
conditions  at  the  next  sample  point,  on  the  basis  of  present  conditions  and  past 
state-vector  aerivitives,  to  introduce  lead  into  the  control  system.  The  control 
generator  is  a  known  (vector)  function  of  the  state  and  independent  variable 
(usually  time).  Its  evaluation  gives  the  optimal  control  for  the  predicted  point. 
The  hold  circuit  supplies  the  optimal  control  over  the  next  time  interval. 
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Development  of  the  control  generator  equation  comprises  the  bulk  of  this  report. 
Section  II  considers  a  two-dimensional  re-entry  trajectory  optimization  problem, 
in  which  total  stagnation  point  heat  (convective  and  radiative)  is  the  function 
minimized.  Terminal  values  of  velocity,  altitude  and  range  are  specified, 
leaving  terminal  flight  path  angle  and  time  unspecified.  After  this  problem  is 
solved,  the  terminal  range  is  extended,  to  obtain  a  design  trajectory,  and  the 
optimal  re-entry  corridor  is  "mapped"  through  the  use  of  a  very  powerful  pre¬ 
dictor  scheme.  The  control  function  for  the  mapped  trajectories  is  used  for 
development  of  the  control  generator  equation.  A  data  generating  program  is 
described  which  supplies  partials  of  the  control  with  respect  to  the  state  vector 
as  well  as  the  control  itself.  This  data  is  punched  on  cards  at  equally  spaced 
values  of  range,  since  range  is  used  as  the  independent  variable  in  the  program. 
It  is  also  shown  in  Section  II  that  the  mapped  trajectories  are  minimizing  tra¬ 
jectories  over  a  large  region  of  solution  space. 

In  Section  III,  the  control  data  is  fitted  as  a  multidimensional  polynomial,  and 
the  fit  is  evaluated,  so  far  as  errors  are  concerned. 


Mechanization  and  simulation  of  the  control  equation  in  the  scheme  of  Figure  1-1 
is  considered  in  Section  IV.  It  is  found  that  the  scheme  produces  ea3onable 
re-entry  trajectories  with  modest  increases  in  total  heat  from  the  optimal 
values  for  most  of  the  trajectories  simulated.  A  few  trajectories  fail  because  of 
control  inaccuracies  in  the  area  of  the  first  dip  into  the  atmosphere,  so  the  region 
of  initial  conditions  for  application  of  the  scheme  must  be  suitably  limited.  The 
region  of  application  is  still  very  large.  The  terminal  condition  errors  for  the 
successful  trajectories  are  found  to  be  largo  (although  tolerable)  and  are  usually 
biased  in  sign  and  magnitude.  It  is  believed  that  these  errors  can  be  greatly 
reduced,  although  no  effort  was  expended  in  this  direction. 

Some  additional  topics  are  considered  in  Section  II.  An  inequality  constraint 
on  the  sensed  acceleration  is  included  in  the  re-entry  trajectory  optimization 
problem.  Several  optimization  methods  arp  used  in  an  unsuccessful  attempt  to 
obtain  a  10-g  optimal  trajectory.  The  most  powerful  method  is  the  predictor 
scheme,  and  this  fails  because  of  a  singular  point  on  the  constrained  subarc. 

A  means  of  isolating  this  singular  point  was  devised,  but  time  limitations  pre¬ 
vented  exploitation  of  the  idea  on  the  computer.  The  optimization  methods  are 
also  extended  to  the  bounded  state  coordinate  problem,  and  the  bounded  brachis¬ 
tochrone  problem  is  solved  on  the  computer. 

Conclusions  and  recommendations  are  presented  in  Section  V. 
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SECTION  II 

OPTIMAL  TRAJECTORY  COMPUTATIONS 


A.  THE  RE-ENTRY  TRAJECTORY  OPTIMIZATION  PROBLEM 


1.  Statement  of  the  Problem 


A  re-entry  path  which  minimizes  the  total  stagnation -point  heating 


of  a  blunt-nose  re-entry  vehicle  is  to  be  found.  The  heating  rate  q  is  the 
sum  of  convective  and  radiative  components 


q  =  qc  +  qr  • 


where 


q  -  cv 
c 


3  n/P 


I  \3/2  .  .12.5 

qr  -  7.5n/_L\  I, - 1 

1  eD  I  lio.ooo/ 

The  state  vector  components  are  the  velocity  v,  flight  path  angle  y, 
dimensionless  altitude  ?  *  -g-,  and  great  circle  range  £.  The  density  p  is 
given  by 


o  -  p  e 
H  vo 


and  values  for  the  constants  are  vehicle  nose  radius  N  =  4  feet,  c  =  2  x  10~  , 
earth  radius  R=  20,  903,  520  feet,  exponential  constant  0  =  1/23,  500  ft.  \ 
and  sea  level  density  p  =  0.0023769  slugs/ft.3. 
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The  equations  of  motion  are 


dv  -S  9 

—  •  - pv  C  <u> 

dt  2m 


dv  S 

-  =  - pv  C.  (u) 

dt  2m  L 


80  sin  y 

7T7T7~ 

v  cos  >  gQ  cos  y 
R(1  +  ?  )  v(l  +  5)2 


d? 

dT 


~R  Sin  7 


(2.6; 


dj  v 

ir =  Tmr ccs  y  • 

S  2  - 1 

The  vehicle  frontal  area  to  mass  ratio  is  — =  0.  5  ft.  -  slug  ,  and  the  lift 

in 

and  drag  coefficients  are  given  in  terms  of  the  control  function  u  as 


CD 

=  CDO  +  CDLcosu 

(2.  7) 

CL, 

=  CLQ  sin  u  , 

(2.  8) 

where 

CDO 

0.  88 

CDL. 

0.  52 

CLO 

-0. 505 . 

Inequality  constraints  include  a  bound  on  the  control  function 

2 

2 

(2.  9) 

U1 

-  u  *  0, 

with  Uj  a  constant,  and  a  pilot's  acceleration  constraint 

B  - 

a  *  0. 

P 

(2.  10) 

G 


In  Equation  (2.  10),  B  is  a  specified  constant,  and  the  pilot's  acceleration  is 
given  by 


where  is  sea-level  gravity  which  normalizes  the  units  of  to  g's. 
Equation  (2.  9)  was  found  necessary  to  produce  initial  trajectories  which 
neither  skipped  out  of  the  atmosphere  nor  dived  in  too  deeply.  It  is  sequen  • 
tially  relaxed  during  the  optimization  process. 


Initial  conditions  are  taken  as  vq  -  35,000  ft/sec,  V q  =  -5.  75  degrees, 

initial  altitude  h  =  400,  000  ft,  and  C  =  zero;  the  terminal  surface  equa- 

o  o 

tions  are 


v(T)  -  Xj  =  0 

h(T)  -  X2  =  0  (2.  12) 

C(T)  -  X3  »  0 

with  the  constants  X^  -  1650  ft/sec,  X 2  =  75,  530  ft,  and  X^  =  979  statute 
miles.  Note  that  the  final  flight  path  angle  and  terminal  time  are  left  un¬ 
specified. 


A  fairlv  detailed  examination  of  this  problem  is  contained  in  Reference  i; 
consequently,  only  a  brief  summary  of  the  results  is  presented  here. 


The  Hamiltonian  may  be  written 

H j  =  q  +  p'f  +  H-jjuj  -  u2]  +  *  ap)  (2.  13) 

where  p  is  the  four -dimensional  multiplier  vector  and  f  represents  the 
right-hand  side  of  the  system  (2.6).  The  Euler-Lagrange  equations,  where 
zero  terms  have  been  omitted,  then  read 
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p4  ■  », |p4  ■  p4  | 


S,1  M2 
0  «  p,  — —  +  P9  — 

1  du  z  du 


da 

2p.«  -  H-o  — 1 B_ 
1  ^  du 


(2.15) 


2.  The  Unconstrained  Subarc 


Both  the  multipliers  Pj  and  are  zero  here,  so  Equation  (2,  15)  is  used  to 
determine  the  control  function.  After  the  substitutions  have  been  made,  the 
resulting  formula  is 


tan  u 


~CLOp2 

_CDLP1V 


(2. 16) 


and  u  is  centered  about  zero  by  the  constraint  (2.  9);  i.  e. , 


-Uj  su  SUj  . 


(2.17) 
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The  minimum -principle  equation  is 


-pjVCpk  cos  u  +  P2cLo  sin  u  *  ’plvCDL  COS  ^  +  P2<'LO  ain  U  *  <2* 

in  which  U  is  any  admissible  value  in  the  range  (2.  17).  The  left-hand  aide 
of  (2.  18)  may  be  considered  as  a  dot  product,  and  the  choice  of  a  unit  vector 
(cos  u,  sin  u)  which  has  minimum  dot  product  with  the  vector  p2^LO^  *8 


sin  u  = 


~CLOP2 


i^~[ 


CDLP1V 


cos  u  = 


CDLP1V 


(2.  19) 


^fLOP^ 


CDLP1V 


This  is  parallel  but  in  the  opposite  direction.  Then,  from  the  signs  of  and 
p2>  assuming  C^q  negative,  it  follows  that: 

If  p  =  0  and  p  >  0,  then  u  -  0 

b  1 


P,  >  0 


Pj  >  0 


0  <  U  <  y 


p2>0 


Pl  =  0 


TT 

u  = 


r>  >0 
K2  w 


p2=  0 


«  s  c\ 

pi  ~ 0 


p  <  c 

1 1 


(2.  20) 


u  =  ±  n  (bang  condition  if  Uj  =  n) 


p2  <  0 


p.  >  0 


-  <  u  <  0 


Po  <  0 


P,  -  0 


n 

u  *  2 


p2<c 


p.  <0 


n  < u  <  -- 
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There  are  no  singular  points  if  and  p2  are  never  simultaneously  zero. 
The  subarc  ends  either  when  (2.  10)  or  (2.  9)  becomes  zero,  or  when  the 
stopping  condition,  the  first  of  (2.  12),  is  satisfied. 


3.  The  Constrained  Subarc  u  =  ±  u. 

Let  0  be  the  angle  defined  by  Equation  (2.  16)  and  the  sign  conventions  given 
by  (2.  20).  Then  substitution  into  the  minimum  principle  equation,  (2.  18),  gives 


coj  (0-  u)  a  cos  (0-  U)  ,  (2.  21) 

which  is  satisfied  if  u  and  0  <  ±  tt,  have  the  same  sign.  The  condition  0  =  tt 
indicates  a  bang.  Furthermore,  substitution  into  (2.  15)  (with  p2  =  0)-  >nd 
some  rearrangement,  gives 


(CLOp2)  +  (CDLP1V) 


2  sin(0-u) 


(2.  22) 


Since  u  and  sin(0-u)  have  the  same  sign,  p,  £  0,  as  required. 

There  are  no  singular  points,  and  the  control  function  is  continuous  at  the 
junction  between  constrained  and  unconstrained  subarcs.  Thus  p^.  from 
(2.  22),  must  start  and  end  with  value  zero,  since  at  such  poinis  u  =  ®. 

Then  the  terminal  surface  is  either  Pj  =  0,  provided  Pj  ^  0,  or  the  stopping 
condition. 


The  constrained  subarc,  over  which  a^  =  B,  is  considered  in  a  separate 
subsection. 
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4.  The  Optimization  Scheme  and  Computer  Results 


It  is  shown  in  Reference  1  that  the  solution  of  systems  of  equations  such  as 
(2.0)  and  (2.  14),  it.  which  u  satisfies  (2.  1 0)  and  (2.  20)  (with  p  j  =  0 ),  or  in 
winch  u  and  p  ^  safisfy  (2.  21)  and  (2.  22),  may  be  written  in  the  form 

X  --  xit,xo.po) 

(2.23) 

P  =-  P<‘.VPo>' 

where  x  and  p  are  the  initial  values  of  the  stale  and  multiplier  vectors 
c  o 

respectively.  It  is  also  shown  in  Reference  1  that  the  necessary  conditions 
which  augment  the  terminal  equations  (2.  12)  are 

P2(T)  =  c 

(2.24) 

H  =  0, 

where  H  is  the  Hamiltonian  (2.  13)  with  the  inequality  constraint  terms  removed. 
When  the  funct.onal  forms  of  the  solutions  (2.  23)  are  taken  into  consideration 
(and  noting  that  xq  is  fixed  for  the  optimization  problem),  the  set  (2.  12)  and 
(2.  24)  become 


v(T.po)  -  Xj  -  0 

?«T-  -V  '  X2  ??  =  0 

C(T,po)  -  X3  =  0 

p^T.Po)  =  0 

H(po)  -  0 


(2.  25) 
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These  five  equations  in  the  five  unknown  quantities  (T,Pq)  may  be  solved  via 
the  Newton-Raphson  method  if  partial  dert.^Uwes  can  be  found.  The  partials 
with  respect  to  T  are  the  time  derivatives  <2.  6)  ana  (2.  14)  evaluated  at  t=  T. 
The  other  partials  are  found  by  integrating  the  system  of  equations 


1  (1*1. 
dt  \3a  j 

d  /  3p\ 
dt  ^3  a  j 


(2. 


in  which  "a"  represents  a  particular  initial  condition.  The  partial  derivatives 
of  the  Hamiltonian  Hj(x,p)  may  change  from  subarc  to  subarc,  and  for  each 
subarc  it  is  known  that 

32H1  32Hj 
3x2  '  dT2 

are  symmetric  matrices  and  that 


>2h, 

“K - > - 

opox 


/32H, 


\~,p/ 


Since  (2.  2(< >  is  a  set  of  2n  linear  first -order,  homogeneous  differential  equa¬ 
tions,  there  is  a  maximal  set  of  2n  independent  (column)  vector  solutions. 
When  initial  conditions  are  taken  as  the  (2n  x  2n)  identity  matrix,  the  first  n 
(column)  vector  solutions  represent  partial  derivatives  with  respect  to  the 
vector  xo>  and  the  last  n  solutions  are  partials  with  respect  to  The  solu¬ 
tions  are  continuous  in  time,  except  possibly  at  corner  points,  where  the 
discontinuities  are  well  defined. 


1 2 


In  the  present  ease,  discontinuities  occur  at  points  t^,  where  the  angle  $ 
becomes  ±  n  ,  or  equivantly,  where  p2  =  0  with  <  0.  Let  rj and 


CijW.  i,  r- 1. 


4,  be  the  elements  of  the 


and 

^  o 


fern 


solution  matrices  respectively.  At  t=  t1 ,  only  the  second  of  Equations  (2.6) 
is  discontinuous  when  u  changes  signs,  and  all  of  Equations  (2.  14)  are 


continuous  (since  p„(t  )  =  01.  According  to  Reference  1,  this  means  that 
i  1  dx 


only  the  second  row  of  -  is  discontinuous  at  t  =  tj,  with 


+  -  2y(t  ) 

tl?1(t.)  =  h,.(t.)  -  - - i~C9.(t,),  j=  1 . 4, 

Zj  1  23  1  p,(t.)  23  1 


(2 


2ivr 


and  (-),  (+)  signifies  values  from  the  left  and  from  the  right,  respectively. 

Now  assume  that  a  path  and  partial  derivative  solutions  have  been  found. 
Then  the  modified  Newton-R;-phson  equations  for  the  system  (2.  25)  are 


p  -1 

- 

— 

-1 

r 

dT 

v(T) 

h.,(T) 

Th2(T) 

r,i,(T) 

VT) 

0 

dp 

Ao 

?  (T) 

'll™' 

T132(T) 

”33,T) 

r,34<T) 

?(T)-X2/R 

_ 

-  -  L. 

C  (T) 

VT> 

M  /rn\ 

n42'T) 

•o  i'T>\ 

’43'  ‘ 7 

n44U/ 

d,\ 

p2(t) 

C2i<t) 

P 

OJ 

VJ 

^23^T' 

C24(T) 

p2<t> 

dP4„ 

0 

fj(0} 

f2(0) 

f3(0) 

f4(0) 

H 

l— 

- 

— 1 

The  zero  in  the  right  hand  vector  corresponds  to  the  stopping  condition,  the 
first  of  Equations  (2.  IT),  which  is  satisfied  by  every  trajectory.  The  Hamil¬ 
tonian  iB  a  constant  for  each  path,  uo  the  last  row  of  the  matrix  contains  its 
partiala  evaluated  at  t  -  0. 


An  optimal  trajectory  for  =  16  degrees  was  obtained  and  is  displayed  in 
Reference  1.  It  appears  here  as  the  first  of  a  family  of  optimal  trajectories 
in  Figures  2-1  through  2-8.  The  second  member,  for  -  20  degrees  (not 
shown),  was  generated  by  using  the  modified  Newton-Raphson  method  with 
the  16-degree  optimal  values  for  (T,  pQ)  as  starting  conditions.  All  other 
members  of  the  family,  obtained  in  a  similar  fashion,  are  displayed  in 
Figures  2-1  through  2-8,  and  values  for  the  optima  criteria  are  given  in 
Table  2-1- 

Figure  2-1  shows  that  the  trajectories  dive  deeper  into  the  atmosphere  and 
skip  higher  as  the  constraint  is  removed.  All  these  curves  end  at  the  same 
point,  as  required  by  the  terminal  surface  equations  (2.  12).  The  first  dive 
produces  higher  pilot's  acceleration  peaks  (Figure  2-2)  but  reduces  the 
secondary  peak.  The  unconstrained  maximum  value  is  20.  5  g‘s.  The  flight 
time  increases  (Figure  2-3)  which  is  a  consequence  of  the  lengthening  skip. 

This  is  apparent  in  the  velocity  curves  (Figure  2-4)  which  tend  to  level  out 
over  the  skipping  portions  of  the  trajectories.  The  flight  path  angles  (Figure 
2-5)  also  show  the  deeper  dive  and  higher  skip.  All  these  curves  apparently 
pass  through  a  common  point,  corresponding  roughly  with  the  bottom  of  the 
first  dip  (see  Figure  2-1).  The  convective  and  radiative  heating  rates  are 
dispiayed  in  Figure  2-6.  They  peak  higher,  and  fall  off  faster,  as  the  con¬ 
straint  is  relaxed.  The  total  heating  rates  of  Figure  2-7  have  the  same  char¬ 
acteristics,  and  show  that  even  though  the  peaks  are  higher,  the  enclosed  area 
becomes  smaller. 

The  optimal  control  functions  are  displayed  in  Figure  2-8,  When  the  16-degree 
trajectory  (which  seems  to  be  in  a  category  of  its  own)  is  excluded,  the  control 
curves  tend  nicely  to  the  unconstrained  trajectory  curve.  They  all  have  a 
"bang"  which  goes  towards  the  endpoint  as  the  constraint  is  relaxed  and,  in 
the  limit,  produces  the  -180-degree  value  of  the  control  function  (the  angle  0 
goes  to  -180  degrees  at  the  bang).  The  first  portions  of  these  curves  show 
that  the  trajectories  are  forced  into  the  atmosphere,  since  positive  control 
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Figure  2-4.  Velocity  Versus  Hange  for  Several  Optimal  Trajectories, 
Parameter  u  ^ 
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Figure  2-6.  Convective  and  Radiative  Heating  RaS.es  Versus  Range 
for  Several  Optimal  Trajectories,  Parameter  Uj 
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Tabic  2-1.  Control  Constraint  uj  ami 
Optimal  Criterion  J 


U1 

(degrees) 

J 

(BTU/ft2) 

U  - 

1G 

27, 334 

25 

2f>,  524 

3  5 

26, 246 

•15 

26, 07] 

55 

25. 951 

180 

25,  73fi 

corresponds  to  negative  lift.  Small  values  of  the  control  also  correspond 
to  maximum  drag,  so  maximum  energy  is  dissipateu.  Before  the  bottom 
of  the  dive,  the  control  functions  all  pass  through  zero  and  then  on  to  the 
maximum  lift  condition  (-90  degrees  for  the  180-degree  optimal).  The 
positive  lift  is  required  for  ranging  purposes  and  is  maintained  for  the  re 
maindc-r  of  the  re-entry  process. 


B.  RE-ENTRY  CORRIDOR  MAPPING  COMPUTATIONS 

It  was  originally  planned  that  optimal  trajectories  having  maximum  sensed 
acceleration  loads  of  10  g's  would  be  used  for  the  nonlinear  optimal  feedback 
control  scheme  of  Section  IV.  However,  a  singular  point  caused  computa¬ 
tional  difficulties.  (The  details  of  the  problem  and  the  method  devised  to 
circumvent  the  difficulty  are  presented  in  another  subsection.)  It  was 
accordingly  decided  to  proceed  using  unconstrained  trajectories  as  a  basis 
for  demonstrating  feasibility  for  the  nonlinear  optimal  feedback  control 
scheme. 
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The  first  step  in  the  corridor  mapping  process  was  to  round  out  the  terminal 
conditions  following  Equations  (2.  12}  to  X2  =  75,000  feet  and  X,  »  1000  miles. 
The  original  values  (75,  530  feet  and  979  miles,  respectively)  were  made 
necessary  by  the  initial  conditions  and  the  almost  ballistic  constraint 
Uj  =  16  degrees.  The  resulting  trajectory  is  displayed  in  Figure  2-9  as  the 
first  member  of  a  family  of  optimal  trajectories  for  which  the  terminal  range 
is  the  parameter. 


1.  Range  Extension 


It  was  decided  to  explore  the  ranging  possibilities  of  the  vehicle  as  the  next 
step  in  the  optimal  re-entry  corridor  mapping  process.  The  predictor  scheme 
presented  in  Appendix  A  greatly  facilitated  the  computations. 

The  differential  equations  are  (2.  6)  and  (2.  14)  (with  p  and  zero),  and  the 
control  u  satisfies  (2.  16)  and  (2.  20),  The  boundary  conditions  are  given  by 
(2.  25)  where  terminal  range  is  the  mapping  parameter.  Differential 
equations  (2.26)  were  integrated  to  obtain  the  partial  derivative  solutions 
for  the  modified  Newton-Raphson  equations  (2.  28).  The  stopping  condition 
for  the  integracions  was  the  attainment  of  a  desired  terminal  time  T  (updated 
after  each  iteration),  sc  the  zero  element  of  the  last  vector  in  (2.  28)  was 
replaced  by  the  term  (v(T)  -  X^).  The  change  in  the  stopping  condition  was 
made  because  it  is  somewhat  easier  and  faster  to  stop  at  a  given  value  of  T 
than  it  is  to  interpolate  for  v(T)  -  Xj.  The  stopping  condition  for  iterations 
was  that  the  maximum  ratio  of  variable-change  to  variable  be  less  than  a 
specified  constant. 

Starting  trajectory  initial  conditions  are  vq  =  35,  000  ft/sec,  >q=  -5.  75  degrees, 
ho  =  400’  000  CQ  =  and  terminal  conditions  are  Xj  =  1650  ft/sec, 

Xg  =  75,  000  ft,  and  Xg  =  1000  miles.  The  next  three  trajectories  were 
obtained  using  the  optimal  Newton-Raphson  method,  described  in  Appendix  B, 
with  terminal  range  increments  of  10  miles.  Succeeding  members  of  the 
family  were  obtained  using  the  Adams -Moulton  predictor  equation 
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(2.29) 


x 

m-r  i 
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h 
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55x'  -  59x'  .  +  37x'  ~  -  9x'  , 

m  m-1  m  -  2  m-o 


in  which  x  is  the  vector  (T,po),  x'  is  the  derivative  of  (T,po>  with  respect 
to  terminal  range,  and  h  is  the  range  increment.  It  is  shown  in  Appendix  A 
that  x  '  is  the  third  column  vector  of  the  inverse  Newton-Raphson  matrix  in 
Equation  (2.  28),  when  (T,  po>  satisfies  Equations  (2.  25). 


The  range  was  extended  to  2020  miles  by  this  process.  Most  of  the  inter¬ 
mediate  predicted  values  of  (T,  pQ)  produced  optimal  trajectories,  which 
shows  the  power  of  the  predictor  scheme.  As  the  upper  limit  on  range  was 
approached,  prediction  gradually  worsened,  anu  it  is  doubtful  that  range 
could  be  extended  much  further  for  the  vehicle  and  initial  conditions  con¬ 
sidered. 


Five  of  the  family  of  trajectories  are  plotted  in  Figures  2-9  through  2-15. 
Figure  2-9  shows  that  the  first  dive  into  the  atmosphere  becomes  shallower, 
as  range  is  extended,  and  that  the  skip  which  follows  becomes  higher  and 
longer.  This  behavior  is  caused  by  the  control  function  (Figure  2-10)  which 
leaves  the  negative  lift  region  (u  >  0),  and  goes  to  the  maximum  lift  condition 
(u  -  -90  degrees)  earlier  in  the  flight  as  range  is  extended.  Less  energy  is 
lost  on  the  first  dive,  as  can  be  seen  in  the  velocity  curves  of  Figure  2-11, 
and  the  decrease  in  the  first  acceleration  peak  of  Figure  2-12.  The  flight- 
path  angle  excursions.  Figure  2-13,  also  become  smaller.  Toward  the  end 
of  the  skip,  the  control  approaches  the  negative  lift  condition,  and,  for  the 
longer  ranges,  produces  negative  lift  to  more  quickly  terminate  the  skip. 

The  secondary  acceleration  peak  rises  with  increasing  range  in  order  to 
dissipate  the  increased  remaining  energy.  There  is  a  minor  sashay  in  the 
paths  near  the  endpoints,  due  to  the  control  passing  through  maximum  lift 
and  (L/D)  on  its  way  to  -180  degrees. 

Figure  2-14  shows  that  total  heat  increases  with  terminal  range,  as  might 
be  expected.  However,  it  is  somewhat  surprising  that  the  total  flight  time 
of  Figure  2-15  first  increases  with  range,  and  then  decreases  for  longer 
terminal  ranges. 
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Figure  2-12.  Sensed  Acceleration  Versus  Range  for  Unconstrained 
Optimal  Trajectories,  Parameter  C  (T) 


timal  Trajectories 


Terminal  Time  Versus  Terminal  Range 
Optimal  Trajectories 


On  the  basis  of  these  results,  it  was  decided  that  the  1500-miie  trajectory 
would  be  used  as  the"  nominal"  trajectory  for  development  of  the  nonlinear 
optimal  feedback  control  scheme.  This  trajectory  represents  a  j  oascsvuble 
compromise  between  total  heat,  peak  acceleration  and  total  length  of  skip. 


2.  Corridor  Mapping  Program 

The  corridor  mapping  problem  considered  is  that  of  sweeping  oi>*  a  reasonable 
region  of  initial  conditions  about  the  "nominal"  trajectory,  and  thereby  ob •• 
taining  a  set  of  optimal  trajectories  covering  the  expected  re-entry  corridor. 
The  terminal  values  of  velocity,  altitude  and  range  are  the  same  for  all  these 
trajectories,  cs  is  the  terminal  value  of  multiplier  and  the  Hamiltonian 
(both  zero),  so  that  the  trajectories  all  belong  to  (he  same  field  of  extremals 
(see  Section  IILA). 

It  is  shown  in  subsection  C  that  time  is  unimportant,  so  far  as  this  problem 
is  concerned,  and  furthermore,  that  all  re-entry  trajectories  maybe  started 
at  a  point  where  initial  range  is  zero.  Hence,  the  only  initial  conditions  which 
need  be  varied  are  vqI  and  5^, 

It  is  convenient,  for  the  mapping  process,  to  integrate  the  system  of  differ¬ 
ential  equations  from  the  terminal  point  to  the  initial  point.  The  predictor 
scheme  may  then  be  used  to  move  the  initial  conditions  over  their  range  of 
variation.  Accordingly,  define  "back  time"  by 


s  -  T-t  . 


(2.30) 


and  note  that 


dx  _  dx  dt 
ds  "  dt  ds 


(2.31) 
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bo  that  the  change  of  independent  variable  changes  only  the  sign  of  the  right- 
hand  sides  of  the  differential  equations.  The  system  of  equations  to  be  integrated 
thus  includes  (2.6).  (2.14)  --  with  p  j  =  p2  “  0  '*  (2.26).  with  opposite 

signs  for  the  right-hand  sides,  and  the  control  u  satisfies  (2.  16)  and  (2.20). 

The  "initial"  conditions  for  the  transformed  equations  '.2.  6)  and  (2.  14)  are 
s=  0  and 


VTXi  =  0 


?T  -X2/R  =  0 


-  X„  -  0 


Pc  ' 

“T 


(2.32) 


withXj  =  1650  ft/s ec,  X2  =  75,000  ft,  and  *  1500  miles.  At  s=  T,  the 
functional  dependence  of  the  solutions  on  the  missing  initial  conditions  may 
be  written 


X=X  (T*  VP3T’P4T) 
p  =  p  (T-  P1T'  *1”  p4T) 


(2.  33) 


so  the  set  of  boundary  conditions  to  be  satisfied  by  each  optimal  trajectory  is 


v(T-PlT‘>T-p3T'p4T 
>(T-plT->T-p3T’p4T 

?(T,PlT'yT'P3T'P4T)  "  X30 
C(T-plT'yT’p3T’P4 
H(PV>T’P3T*P4T)  =  0 


X10  *  0 


X20  -  0 


3C/R  =  0 

=  0 


(2.34) 
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XJ0>  X20  and  X^^  are,  of  course,  the  initial  conditions  at  t«  0  (or  terminal 
conditions  at  ss  T)  to  be  varied  in  the  mapping  process. 


Let  n_(s)  and  £^(8),  i.  j=  1,  ....  4,  be  the  solutions  of  the  transformed  system 
of  Equations  (2.  26)  with  initial  conditions 


~0  0  0  o' 

I 

10  0  0 

0  10  0 

:jo)\ « 

0  0  0  0 

0  0  0  0 

L  11  j 

0  0  10 

0  0  0  0 

0  0  0  1 

(2.  35) 


Then  the  Newton-Raphson  equations  for  system  (2.  34)  are 


r*  ~ 

dT 

v(T) 

VT) 

^12<T> 

"l3<T> 

h14(T) 

-1 

dp, 

T 

y(T) 

t,21(T> 

r,22*T* 

vT> 

h24(T) 

y  (T)-X20 

dy-p 

-  - 

?  (T) 

\(T) 

t132<T) 

VT) 

T134(T> 

?(T)-x30/r 

dp, 

T 

C<T) 

VT> 

VT> 

TWT> 

C(T) 

dp4 

T 

0 

fj<0> 

-P2(0) 

f3(0) 

14<C> 

H 

(2.36) 


Tiie  first  three  columns  of  the  inverse  matrix  may  be  identified  as  the  deriva¬ 
tives  of  (T,  piT>  yT,  P3t,  p4„J  with  respect  toXj0.  X2£).  and  X3Q/R,  res¬ 
pectively.  for  the  predictor  equation  (2.  29). 

Twenty-six  trajectories  spanning  the  re-entry  corridor  were  generated  using 
the  optimal  Newton-Raphson  method  and  the  predictor  scheme.  The  initial 
conditions  and  total  heat  for  these  and  the  "nominal”  trajectory  (at  t=  0)  are 
summarized  in  Table  2-2.  Trajectory  2  was  obtained  first.  The  first  three 
trajectories  were  determined  by  the  optimal  Newton-Raphson  method,  using 
&X,q  =  +50  ft/ sec  as  the  increment.  The  predictor  scheme  then  generated 
the  remaining  trajectories.  The  predictor  scheme  was  used  exclusively  to 
obtain  trajectory  3,  since  back  derivatives  were  available  from  the  trajectories 
leading  up  to  trajectory  2. 
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Table  2-2.  Initial  Conditions  and  Total  Heat  for  27  Optimal 
Trajectories  Obtained  During  the  Corridor 
Mapping  Process 


Initial  Sign  of  Change 

Initial  Flight  Path  Initial  From  Trajectory  1 

Trajectory  Velocity  Angle  Altitude  AX  |  AX„n  I  AX,„ 
Number  (ft/sec)  idegrees)  (It)  10  20]  30 


35000 
35750 
34250 
35000 
3  5000 
35000 
35000 
35750 
35750 


35750 
34250 
34250 
3  50C0 
35C00 
35000 


35000 

35600 

35600 

35600 

35600 

34400 

34400 

34400 

34400 


-5.  75 

400, 00C 

-5.  75 

400, 000 

-5.  75 

400,  000 

-5.  30 

400, 000 

-7.  65 

400,  000 

-5.  75 

440, 000 

-  5.  7  5 

360,  000 

-5.  30 

400, 000 

-7.  65 

400, 000 

-5.  30 

400, 000 

-7.  65 

400, 000 

-5.  75 

440, 000 

•  5.  75 

360, 000 

-5.  75 

^40, 000 

-5.  75 

360, 000 

-5.  30 

427, 50u 

-5.  30 

i 

360, 000 

-7.  65 

440, 000 

-7.  65 

360, 000 

-5.  40 

430, 000 

-5.40 

370, 000 

-7.  35 

430, 000 

-7.  35 

370, 000 

-5.  40 

430, 000 

-5.  40 

370, 000 

-7.  35 

430, uuO 

-7.35 

370, 000 

Total 

Heat 


28. 872 
31, 888 
26. 252 
29, 278 
28, 66G 
29, 032 
28, 757 
32. 280 
32, 116 


26, 657 
25, 792 
32. 135 


29.  017 


32, 

444 

3 1 

374 

31, 

110 

31, 

375 

27, 

326 

26,  216 
26, 368 


Trajectories  4  and  5  were  generated  in  a  similar  fashion,  using  =  0.  025 

degree.  It  was  found,  however,  that  as  the  flight  path  angle  became  steeper 
(going  toward  trajectory  5),  the  increment  could  he  increased  to  =  0.  1 

degree  with  very  little  degradation  in  the  performance  of  the  predictor  scheme. 

Trajectories  6  and  7  were  obtained  similarly,  and  it  was  determined  that 

=  25C0  ft  was  a  satisfactory  increment  for  the  altitude  mapping  process. 

From  convergence  characteristics  of  the  optimal  Newton-Raphson  process,  it 
was  found  that  velocity  changes  were  easiest  to  obtain,  and  that  flight  path 
angle  changes  wore  hardest  to  satisfy.  Consequently,  trajectories  8  through 
19  were  generated  using  the  most  easily  changed  single  parameter.  For 
example,  the  trajectory  for  which  '‘20  =  ‘5‘  30  degrees  was  the  starting  path 
for  the  generation  of  trajectories  8  and  10,  and  the  mapping  process  look  place 
over  Xjp. 

Trajectories  20  through  27  were  generated  in  a  similar  fashion,  except  that  inter¬ 
mediate  maps  over  X^0  were  performed  first,  followed  by  maps  over  Xjq.  Thus 
for  example,  the  trajectory  for  which  X^q  =  -5.  40  degrees  was  the  starting  tra¬ 
jectory  to  generate  paths  for  which  =  370,000  feet  and  430,000  feet  respec¬ 
tively.  The  Xj0  map  for  these  paths  produced  trajectories  21,  25,  20  and  24. 

The  ease  of  attainment  of  the  various  paths  varied  considerably  as  a  function 
of  parameter  values.  As  pointed  out  above,  the  step  size  for  the  XgQ-map 
could  be  increased  to  =  0.  1  degree  as  the  flight  path  angle  became 

steeper.  When  going  in  the  opposite  direction,  however,  the  smaller  step 
size  was  required.  Trajectory  17  was  generated  with  some  difficulty,  al¬ 
though  the  process  became  easier  as  initial  altitude  decreased.  In  trying  to 
increase  altitude  to  440,000  feet  for  trajectory  16,  the  predictor  scheme 
failed  to  go  beyond  X^q  =  427,  500  feet. 
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In  most  cases,  the  velocity-map  was  easily  accomplished.  Optimal  Newton- 
Raphson  many  times  converged  in  one  step,  and  the  predictor  scheme  (once 
started)  almost  always  predicted  optimal  trajectories.  Some  difficulties 
were  experienced  in  trying  to  obtain  trajectories  20  and  24  (the  worst  cases 
were  always  shallow  initial  flight  path  angle  and  increased  altitude,  as  might 
be  expected).  A  trajectory  for  XjQ  -  34,950  feet  was  obtained  after  several 
optimal  Newton-Raphson  iterations,  and  further  such  steps  would  have  been 
wasteful  of  computer  time.  A  two  point  predictor  formula  from  Reference  2 
was  used  to  predict  trajectories  above  and  below  the  two  members  of  the 
\elocity-map  family  of  paths.  With  these  improved  estimates,  the  optimal 
Newton-Raphson  scheme  converged  rapidly,  and,  thereafter,  the  predictor 
scheme  predicted  optima)  trajectories  at  each  step.  This  again  shows  the 
power  of  the  predictor  scheme  and  indicates  that  lower -ordered  predictor 
formulas  might  well  be  of  value  during  the  mapping  process. 

In  all  cases,  the  optimal  Newton-Raphson  scheme  was  quite  sensitive  to  the 
weighting  factors  used  to  multiply  Equations  (2.34)  (see  Appendix  B).  It  was 
found,  however,  that  a  single  set  of  weights  could  be  used  when  mapping  over 
a  given  variable.  The  numbers  used  are  given  in  Table  2-  3. 

Table  2-3.  Weighting  Factors  for  the  Optimal  Newton-Raphson 
Method  Used  for  the  Mapping  Variables. 


W  oi 

Y 

'MO 

Variable 

X20 

V 

'*30 

W1 

10 

lx  10'2 

i 

W2 

4 

1  x  10 

1 

1 

4 

1  x  10 

W3 

1  x  10 

1  x  10'2 

10 

W4 

1  x  10'3 

1  x  10'3 

W5 

— 

i 

1 

1 

v 


i 

i 

i 


i 

! 


1  ' 
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Trajectories  2  through  7  are  compared  with  the  "nominal"  1500-mile 
trajectory  in  Figures  2-16  through  2-27,  It  is  seen  that  trajectories  2  and 
3  are  quite  symmetric  about  the  nominal  path,  except  for  areas  where  the 
curves  cross  each  other.  Larger  differences  from  nominal  values  are  observed 
for  trajectories  4  and  5  (Figures  2-20  -  2-23),  and  particularly  for  trajectory  5. 
Note  that  the  steeper  initial  flight  path  angle  produces  a  large  skip  and  high 
acceleration  peaks.  Somewhat  similar  results  are  observed  for  trajectories 
6  and  7  (Figures  2-24  -  2-27).  Here  however,  a  lower  initial  altitude  produces 
a  longer  skip  and  higher  acceleration  peaks.  The  control  functions  for  all  27 
optimal  trajectories  are  plotted  in  Figure  2-28.  Note  that  the  characteristics 
are  similar  to  those  of  Figure  2-10.  The  control  functions  for  those  trajec¬ 
tories  displaying  large  skips  assume  their  minimum  values  earlier,  and  then 
rise  higher  than  those  for  the  better  behaved  paths.  A  good  coverage  of  the 
control  region  of  interest  was  obtained  during  the  mapping  process. 


C.  THE  DATA  GENERATING  PROGRAM 

The  data  to  be  fit  consists  of  the  control  function  u  for  the  27  optimal  tra¬ 
jectories,  and  two  multipliers  like  p^  and  p^.  The  multipliers  are  included 
in  hopes  that  an  arctangent  relationship  such  as  (2.  16)  will  produce  an  over¬ 
all  reduction  in  the  size  of  the  control  function  fit.  Additional  data  includes 
the  partial  derivatives  of  u  and  the  two  multipliers  with  respect  to  the  state 
vector. 


1 .  Statement  of  the  Problem 

The  problem  is  transformed  to  an  equivalent,  although  somewhat  reduced 
and  more  convenient,  form  through  a  series  of  transformations.  These 
transformations  are  carried  out  in  Appendix  C  and  the  results  are  presented 
here.  The  new  set  of  differential  equations  corresponding  to  (2.2)  and  sys¬ 
tem  (2.  6)  are: 
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Figure  2-16.  Velocity  Versus  Range  for  Optimal  Trajectories 
1,  2,  and  3 


Versus  Range  for  Optimal  Trajectories 


2  and  3 


033)1330)  VHB» 
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Figure  2-21.  Flight  Path  Angle  VersuB  Range  for  Optimal  Trajectories 
1.  4  and  5 


Figure  2-22.  Altitudf  Versus  Range  for  Optimal  Trajectories  1,  4  and  5 


Figure  2-23.  Sensed  Acceleration  Versus  Range  for  Optimal  Tra¬ 
jectories  1,  4  and  5 


ec 


dq  Uj  (R+h)  p  l/2  V  a,  (R+h)  /  p  \3/2  /  V  \  4 


<os  y  p , 


dV  pV(R+h)  a4  tan  y 

dz  '  cos  >  CD(u)  +  (R+h) 


ac  p(R+h) 


<2.37) 


"Hz  "  a5  +  cosy  CL(u)  +  V(R+h) 


=  a5  (R+h)  tan  y 


dt  a  (R+h) 

—  s  iTT - 

dz  V  cos  y 


with  constants  defined  by 


a!  1  V^ras 

a9  =  7.  5  x  10*4N  a. 


a4  =  2a7 


<2.38) 


Cog  R. 
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The  new  independent  variable  is  range  measured  in  miles,  which  is  related 
to  C through 


c  3  =  5,280 


and  the  new  velocity  and  altitude  are  given  by 

, .  2 
V  =  v 

h  =  R  ?  . 

The  terminal  conditions  corresponding  to  (2.  12)  are 


V(zT)  -  X“  =  0 


h(zT)  -  X2 


ZT  "  ^3^  c3  - 


(2. 39) 


(2.40) 


(2.41) 


The  new  Hamiltonian  is 


h2  =  go  +  p  g- 


(2.42) 


where  gQ  is  the  right-hand  side  of  the  first  of  (2.  37),  g  is  a  four-dimensional 
vector  denoting  the  remaining  right-hand  sides,  ind  P  is  the  new  multiplier 
vector.  The  Euler-Lagrange  equations  are 
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dg 

8°  4- 

dgl 
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"5V  + 

P1 

3V 
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3g2 
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r»  °  , 
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+  P2  ay  + 

P3  ay 

P4 

by 

bg 

°  + 

P 

^1 
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D 
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dz 

3h 

P1 

3  h 
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P4 

ah 

(2. 43) 
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The  control  function  satisfies 


u  =  tan 


,CLOP2 

2CDLP1V 


(2.44) 


with  the  sign  conventions  (2.  20).  Since  P4  is  a  constant  and  terminal  time  is 
an  unspecified  quantity,  the  transversality  conditions  require 


P4  -  0. 


(2.  45) 


Thus,  the  time  equation,  the  last  of  (2.  37),  is  unimportant  in  the  problem, 
and  the  set  of  dependent  variables  is  reduced  to  (V,  y,  h). 

Other  boundary  conditions  could  be  derived.  However,  this  problem  is 
equivalent  to  the  problem  of  subsection  A  (when  ^  K  -  0)  with  equivalence 
relationships  (2.39),  (2.40),  and,  as  shown  in  Appendix  C 

P  -_1l 

1  2v 


(2.46) 


P„  =  -H 


H2  =  '  c  3P4‘ 


2.  The  Partial  Derivatives 

The  system  of  variational  equations  like  (2.26)  for  the  new  system  of 
differential  equations  may  be  written 
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Let  the  fundamental  system  of  solutions  to  (2,47)  be 


tt(z)  = 


dx(z) 

dP(z) 

■^7 


dr.(z) 

o 


n(C)  -  I, 


and  furthermore,  let  the  elements  jf  the  submatrices  be  denoted  by 


bx(z) 

V1’ 

a  x(z) 

Vz) 

*Xo 

,  j=  1,  2,  3,  g  p 

o 

bF(z) 

a  x 

° 

. 

a  p(z) 

,  i,  j-  1,  2,  3,  = 

"  o 

i  *  1,  2,3 
j  =  4,  5,6 

i  «  1,2,3 
3  =  4,5,6 


(2.47) 


(2.48) 


(2.49) 


An  optimal  trajectory  in  the  neighborhood  of  an  optimal  trajectory  will  satisfy 
the  perturbed  equations 


Cx(z) 


6P(z> 


d  x(z) 
Tx 
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fi*o 
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3x(z) 


CP 
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a  p(z) 

bX 

o 


Cx 
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ap<z) 

3Po 


CP 

O 


<2.  50) 


at  arbitrary  values  of  range,  and  linearized  versions  of  the  boundary  c  editions 
(2.41),  (2.45)  and 


at  the  terminal  point.  These  may  be  written 


— 

_  _ 

—  _ 

_  _ 

CV(zT) 

0 

r'n(i.T)  r^12(zT)  V13'ZT* 

ftVo 

nl4(zT)  ri15(zT)  ^  16^ZT^ 

6pio 

6h(z^,) 

= 

0 

= 

r,3l<zT>  ’132(ZT) 

+ 

T134(zT>  r35(zT)  t136(zt) 

6P20 

CP2(zT) 

0 

C2i(z.T)  C22(zT)  C23(zT) 

Ch 

L  j 

C24(zt>  C25(ZT)  C26(zT} 

5P3h 

(2.  52) 


since  equation  (2.  45)  and  the  last  or  (2.  41)  are  al.ay,  satisfied  in  the  reduced 
system.  Let  (2.52)  be  represented  by  - "“v't: — 


A  Cx  +  E  CP  =  0 

o  o  ' 

so  that 


«PC  =  -B_1A  6x  =  K  Ci  . 
u  o  o 


<2.  54) 
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Note  that  B  is  the  Newton-Raphson  matrix  in  the  reduced  system,  which  must 
be  nonsingular.  It  is  shown  in  Appendix  E  that  if  B  is  nonsingular,  the  path  is 
normal  according  to  the  definition  of  normality  in  Reference  7, 

When  (2.  54)  is  substituted  into  (2,  50)  the  result  is 


Under  the  assumption  that  the  first  matrix  of  (2.  55)  is  nonsingular,  it  follows 
that  the  partial  derivatives  of  the  multipliers  with  respect  to  the  state  are 


12.  56) 


The  partial  derivatives  of  the  control  are  derived  from  Equations  (2.44)  and 
(2,  5C)  with  the  result  being 


[~  3  u 

j^dv  ■ 


0, 


(2.  57) 


Several  comments  should  be  made.  The  process  described  above  constructs 
a  particular  field  of  extremals  in  the  neighborhood  of  an  extremal  trajectory 
if  the  inverse  matrix  in  (2.  5G)  exists  at  each  point  along  the  path  (see  Reference 
7,  pp  237-240  and  Section  IlLA  following).  By  construction,  however,  the 
matrix  is  singular  at  the  endpoint.  In  Bpite  of  this  singular  poi^t.  it  can  be 
shown  that  the  results  do  constitute  a  field  when  the  endpoint  ip  excluded,  so 
long  as  the  inve-’se  matrix  in  (2,  56)  exiets  at  oil  other  points  along  the  extre¬ 
mal  path.  The  27  optimal  trajectories  extend  this  field  over  a  finite  region  of 
space.The  partial  derivatives  (2.57)  may  be  UBed  to  construct  a  "neighboring 
optimal"  linear  feedback  control  scheme,  as  in  Reference  8,  about  any  one 
of  the  optimal  trajectories, 
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3 .  Computer  Results 


The  system  of  differential  equations  in  the  computer  program  includes  (2.37) 
and  (2.43),  in  which  the  last  equation  of  each  set  is  omitted,  and  the  control 
function  satisfies  (2.44)  and  (2.20).  The  equivalence  relationships  (2.40) 
and  (2.46)  establish  initial  conditions  from  the  results  of  the  mapping  pro¬ 
gram.  System  (2.47)  is  also  integrated  to  obtain  the  solutions  (2.48).  Since 
the  initial  conditions  for  the  multipliers  may  be  somewhat  in  error,  Newton- 
Raphson  equations  arc  included  to  adjust  them  to  their  proper  values.  When 
the  process  has  converged,  the  K-matrix  of  Equation  (2.  54)  is  computed,  and 
the  iniita!  conditions  for  system  (2.47)  are  changed  to 

’  I  I 

n  t0)  = 

_  K  -I 

for  the  last  pass.  The  first  three  column  solutions  are  then  the  matrices  of 
Equations  (2.  55),  and  the  last  three  solutions  are  used  for  an  auxiliary 
sufficiency  condition  computation  (see  subsection  D).  On  the  last  pass  the 
control  u,  multipliers  Pj  and  P^.and  the  partial  derivatives  [from  (2.  56)  and 
(2,  57) j  are  punched  on  cards  at  specified  values  of  range.  In  the  present 
case  101  sets  of  data  are  obtained  spaced  at  range  increments  of  15  miles. 

At  the  last  point  (z  =  15U0  miles)  the  partial®  are  omitted. 

In  computing  the  partiuls  at  tne  output  points,  the  inverse  matrix  of  Equation 
(2.  56)  is  obtained  by  inverting  the  solution  in  (2,  55),  This  does  not  increase 
computing  time  excessively,  since  the  matrix  Is  small  (3  x  3).  In  larger 
systems,  however,  one  may  wish  to  solve  the  nonhomogeneous  Ricat.ti  equa¬ 
tion  derived  la  Appendix  D.  and  thereby  obtain  the  solutions  (2.  56)  directly. 

The  partial  derivatives  of  the  control  function  for  the  nominal  trajectory  arc 
plotted  in  Figure  2-29.  The  curve  appears  to  be  rather  uninteresting  over 
most  of  the  trajectory,  until  it  )s  remembered  that  V  is  a  very  large  quantity 
over  this  region.  The  hump  in  the  curve  corresponds  approximately  to  the 
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Figure  2-29.  Partial  Derivatives  of  the  Control  With  Respect  to  V, 
y  and  h  as  Functions  of  Range  for  Trajectory  1 


bottom  of  the  first  dive  (see  Figures  2-16  through  2-18),  where  control  changes 
drastically  influence  the  remainder  of  the  path.  Beyond  this  point,  the  partials 
are  small,  corresponding  to  the  skip  in  the  trajectory.  The  partials  all  head 
toward  +  <*>  at  the  end  of  the  trajectory  (by  construction).  Since  the  partials 
can  be  used  as  time  varying  gains  in  a  linear  feedback  control  scheme  (as  in 
References  1  and  8)  one  might  wonder  why  they  are  not  negative,  to  corres¬ 
pond  to  negative  feedback.  The  reason  for  this  is  that  positive  control  pro¬ 
duces  negative  lift.  One  can  intuitively  justify  that  the  partials,  acting  as 
feedback  gains,  produce  the  proper  trajectory  changes,  at  least  toward  the 
end  of  the  path. 

Figures  2-30  and  2-31  are  plots  of  the  multipliers  and  for  trajectories 
1  through  7.  They  are  included  to  show  the  nature  of  the  variations  over  the 
optimal  re-entry  corridor. 

SUFFICIENCY  CONDITIONS 

The  sufficiency  tests  described  in  this  subsection  are  tailored  to  problems 
like  the  re-entry  problem  stated  above  in  which  the  Lagrange  form  of  the 
optimization  problem  is  considered  and  the  terminal  equations  are  quite 
simple.  The  tests  can.  of  course,  be  extended  xo  other  cases,  such  as  prob¬ 
lems  with  more  complex  endpoint  equations,  the  addition  of  a  function  of  the 
endpoints  to  the  function  to  be  minimized,  and  the  inclusion  of  inequality 
constraints  in  the  problem  statement. 

1 .  Sufficiency  Conditions  for  a  Relative  Minimum 

Bliss’  theorem  85.  1  (Reference  7,  page  241)  forms  the  basis  for  the  relative 
minimum  sufficiency  test.  This  theorem  requires 
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Multiplier  P  for  Several  Optimal  Trajectories 


(1)  A  field  in  which  the  extremal  is  imbedded 

(2)  The  strengthened  version  of  the  Weierstrass  condition 

H  (u)  <  H(U)  ( 2 .  59) 

(3)  The  second  variation,  evaluated  with  field  elements,  to  have 
a  proper  minimum  at  the  ends  o.  the  extremal. 

If  these  three  conditions  at  e  satisfied  (for  a  normal  extremal  without  corners) 
then  the  path  is  a  strong  relative  minimizing  path  in  the  tense  of  Bliss' 

Theorem  82,  1  (Reference  7,  page  235). 

For  application  to  the  re-entry  problem,  consider  the  formulation  of  sub¬ 
section  C.  Condition  (2)  is  satisfied  since  the  control  (2.44)  using  sign  con 
ventions  (2.  20)  establishes  the  absolute  minimum  of  the  Hamiltonian  (with 
respect  to  u)  and  no  other  control  U  gives  the  same  value.  The  problem  is 
normal  s  ince  the  Newton-Raphson  matrix  E  of  Equation  (2.  53)  is  nonsingular. 

It  is  shov/n  in  Appendix  E  that  B  nonsingular  implies  normality  in  the  sense 
of  Reference  7. 


A  field  is  constructed  through  the  use  of  lemma  84.  2  of  Reference  7  (pp  238 
240).  This  requires  that  the  path  first  be  nonsingular.  Lemma  87.  1  of 
Reference  7  (page  247)  shows  that  the  path  is  nonsingular  if  the  strengthened 
Clebsch  condition  (Theorem  78.  2  of  Reference  7)  holds.  With  a  single  control 
function  this  requires 


>  0 


(2.  60) 


over  the  path,  which  is  true  for  the  problem  under  consideration.  Then  a 
field  is  constructed  if  a  conjugate  system  of  solutions  U,  V  (both  nxn  matrices) 
to  the  accessory  equations  (2.47)  can  be  found,  with  U  nonsingular  everwhere 
along  the  extremal.  By  definition,  a  conjugate  system  of  equations  satisfies 

u'v  =  V'u.  r  (2.61) 
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A  suitable  choice  for  l!  and  V  is 


U 


-V 


3x  3x 

IT'  '  Tp 

o  o 

dp  3  P 

3x  ’  3P 
o  o 


nll<z)  +  r,i2(z)  [K'Ii 
n21(z)  +  n22(z)  !K'J] 


(2.  62) 


which  accounts  for  the  second  half  of  the  initial  conditions  (2.  58)  for  Equations 
(2.47).  The  nnotation.  adopted  here,  splits  the  fundamental  solution  matrix 
(2.  48)  into  the  four  submatrices  cated  in  (2.62).  One  can  show  that  the 
choice  (2.62)  forms  a  conjugate  system  because  of  the  easily  proved  relation¬ 
ships 
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Then  a  field  is  constructed  if  U  of  (2.62)  is  nonsingular  at  each  point  along 
the  extremal. 


Several  comments  should  be  made.  Although  the  results  (2.  62)  and  (2.  63) 
are  presented  for  the  problem  of  subsection  C,  they  hold  in  general.  The 
minus  sign  in  (2.  62)  arises  because  of  the  use  of  the  minimum  principle 
instead  of  the  maximum  principle,  as  explained  in  Appendix  E.  Since  the 
original  problem  included  u  as  the  derivative  of  an  additional  state  coordinate, 
one  might  wonder  whether  the  field  should  be  constructed  in  (n+2)  dimensions 
rather  than  in  (n+ 1)  dimensions.  It  is  shown  in  Appendix  E  that  the  additional 
dimension  need  not  be  considered.  The  sufficiency  proof  is  complete  at  this 
point  for  fixed  endpoint  problems. 

^Matrix  K.  from  Equation  (2.  54)  is  taken  as  zero  in  (2.  58).  It  is  only  necessary 
to  show  that  there  is  a  conjugate  system  U,  V  with  det  U  /  0. 
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In  ihe  computer  program,  the  de.er tmnant  of  U  was  computed  at  each  output 

point  (z- 0  to  1500  miles  with  increments  izz  -  15  miles).  The  determinant 

started  at  unity  (z  -  0),  and  rose  to  very  large  values  as  z  increased  (about 
29 

10  maximum)  for  each  of  the  27  trajectories  considered.  Because  of  the 
smoothness  of  the  problem  solutions,  it  was  concluded  that  the  determinant 
did  not  vanish  or  become  negative  between  output  points,  and  hence,  that  fields 
were  indeed  constructed  for  each  trajectory. 


Condition  (3),  above,  is  based  upon  lemma  (87.2)  of  Reference  7  (page  247). 
This  gives  generalizea  expressions  for  the  second  variation  test  for  problems 
with  separated  end  conditions.  There  is  no  contribution  to  the  second  variation 
at  the  initial  point,  since  conditions  at  this  point  are  all  fixed.  Hence,  for  the 
type  of  problem  considered,  it  is  required  that 


a  t 

o 


p ( z T )  x('/T) 
U'(zT)  p(zT) 


p  '(zT)  U(zT) 

ao 

U'(zT)  V(zT) 

a 

(2.  64) 


for  all  <aQ,  a)  satisfying 


xi(zT>aQ  +  lh(/  ,)a  -  0,  i-  1 .....  r. 


(2.  65) 


Expression  (2.64)  is  the  second  variation  (as  derived  in  Reference  1)  evaluated 
at  the  terminai  point  with  field  elements,  and  (2,65)  represents  the  linearized 
terminal  conditions  with  subscript  i  on  U  representing  the  ith  row  vector  of  U. 


In  the  re-entry  problem  the  terminal  value  of  range  is  also  specified,  which 
implies  aQ  -  0,  and  accordingly  reduces  the  size  of  expression  (2.64).  The 
remainder  of  the  terminal  conditions  (2.  41)  may  be  written  in  the  linearized 
form 


U(zT>  a 


(2.  66) 
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where  the  middle  row  of  U  is  included  for  the  unspecified  variable  y,  and  A 
is  an  arbitrary  scalar.  When  {2.66)  is  inserted  in  (2.64)  (with  ao  =  0)  there 
results 


-  0  A  oj 


V(zT)U‘1(zT) 


>  0. 


(2.67) 


Thus,  the  center  element  of  the  product  matrix  in  (2.  67)  must  be  negative  to 
satisfy  condition  (3).  It  is  negative  for  each  of  the  27  trajectories  considered. 
Hence,  the  sufficiency  conditions  are  all  satisfied,  and  each  of  the  27  trajectories 
is  a  stong  relative  minimizing  trajectory. 


2.  An  Absolute  Minimum  Test 

The  test  given  here  does  not  establish  global  sufficiency.  It  does,  however, 
allow  a  large  region  of  solution  space  to  be  examined  for  other  solutions  to 
the  optimization  problem.  The  method  is  not  applicable  to  fixed  end-point 
problems.  The  idea  is  simple:  replace  one  of  the  transversality  conditions  by 
a  new  terminal  equation  in  which  a  parameter  is  included.  Obtain  a  set  of 
soultions  to  this  problem,  as  functions  of  the  parameter,  and  examine  the  set 
for  satisfaction  of  the  omitted  transversality  condition. 

To  illustrate  the  method,  consider  the  problem  statement  of  subsection  A,  and 
in  particular,  the  set  of  necessary  condtions  (2.  25).  The  fourth  of  these  is 
replaced  by  the  equation 

y(T,po)  T  =  0  (2.68) 

in  which  1'  is  the  parameter,  and  solutions  as  functions  of  T  are  to  be 
examined  for  satisfaction  of  p  V(T)  =  0.  The  Newton  -Raphs on  equations  for 

the  new  system  of  equations  may  be  written 
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According  to  Appendix  A,  the  derivatives  of  the  initial  condition-i  with  respect 
to  the  parameter  r  are  in  the  second  column  of  the  inverse  matrix.  These 
are  used  with  the  predictor  equation  (2.  29)  to  establish  the  family  of  solutions. 

An  unconstrained  optimal  trajectory  with  terminal  range  of  1450  miles 
was  used  for  the  test.  The  results  certainly  hold  for  the  1500  mile 
trajectory.  During  computer  runs  it  was  found  that  an  increment 

AT  =  1/2  degree  initially  produced  a  member  of  the  family  of  solutions  at 
each  step.  Prediction  gradually  worsened,  however,  as  the  extremes  of  F 
were  approached,  and  it  is  unlikely  that  further  extension  of  F  can  be  accom¬ 
plished. 

The  total  heat  J  for  the  family  of  trajectories  is  plotted  in  Figure  2-32. 
There  is  only  one  minimum  for  the  range  of  F  given,  and  it  is  noted  that  the 
minimum  is  quite  insensitive  to  large  variations  in  terminal  flight  path  angle. 

From  Hamilton-Jacobi  theory,  it  is  known  that-p^(T)  is  the  slope  of  Figure 
2-32.  The  slope  is  zero  for  the  optimal  trajectory,  and  becomes  very  large 
(in  absolute  value)  at  the  extremes  of  the  curve,  so  -p.-,(T)  is  never  again 
zero  over  the  range  of  F.  It  is  concluded  that  the  original  optimal  trajec¬ 
tory  is  the  absolute  minimizing- trajectory ,  at  least  over  the  obtainable  range 
of  the  parameter  F. 
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Tigure  2-32.  Absolute  Minimum  Test  Results  -  Total  Heat  J  Ver 
Terminal  Flight  Path  Angle  T 
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E.  THE  CONSTRAINT  a  =  B  AND  EXPERIENCES  WITH  VARIOUS  OPTIMIZATION 
SCHEMES  p 


1.  Constrained  Subarc  Equations 

Consider  the  problem  of  subsection  A,  and  include  the  constraint  equation 
<2.  10)  in  the  problem  formulation.  When  equality  holds  in  (2.  10),  the  resulting 
equation  is  used  to  determine  the  control  u.  This  is  squared  and  rearranged  to 
obtain 

a  cos^u  +  2b  cos  u  +  c  s  0  (2.  70) 


where 


a  C  ^  _  ('  2 

LDL  CL.O 


CDOCDL 


(2  71) 
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Upper  and  lower  limits  for  application  uri'  found  by  substituting  u  0.  i  v  into 
(i.  70)  and  expanding  in  the  vicinity  of  these  points.  It  is  found  that  the  inquality 


(cix>  f  * 


* *  a 

•» 

Sp  v“ 


*  (CDf)  '  CDL)  >  0 


(2  72) 


v  a 

i"uip  lioid.  The  expression  fur  -y^-  is  proportional  to  am  u,  which  is  aero 
at  u  0,  Iff,  These  points  are  thus  singular  points,  bo  equality  must  be 
excluded  In  (2.  72).  Assuming  this,  Equation  (2.  70)  may  be  solved  to  obtain 


Cos  u  »  — 
a 


1  •  V77^ 


(2.  72) 


where  the  omitted  root  lulls  outside  the  range  |  cos  u  |  4  i.  It  is  easily  shown 
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that  the  term  in  the  square  root  is  positive  by  substituting  the  upper  limit  of 
(2.  72)  into  the  expression  and  evaluating. 

The  minimum  principle  again  takes  the  form  (2.  21),  except  that  this  time 

-  iu  i  <  0<  | u  I  .  .  (2.  74) 

Otherwise,  <t>  and  u  would  be  identical.  7’hen,  once  again,  $  /  0  and  u  must 
have  the  same  sign,  and  <£  =  0  is  the  bang  condition. 

The  multiplier  determined  from  Equation  (2.  15)  (with  p.j  r  0).  When  all 

the  substitutions  arc  made,  it  turns  out  that 


p 


2 


v 


sin(u-«  ) 


sin  u 


(2.  75) 


Since  sin  (u-i/l)  and  sin  u  have  the  same  signs,  p^  x  0  as  required.  Again,  u 
is  continuous  at  junction  points,  so  p^  must  start  and  end  at  zero;  and  p^  =  0 
with  p,,  to,  describes  the  subarc  terminal  surlace,  unless  the  stopping  condi¬ 
tion  is  satisfied  first. 

It  is  convenient  to  include  H  as  an  additional  parameter  in  the  problem,  and  to 
introduce  the  new  terminal  equation 

D  -  B  *  0.  (2.  70) 

o 

H  is  the  final  maximum  g-level  (usually  token  as  10  g's)  in  the  constraint 
equation  (2.  10),  and  the  solutions  now  depend  upon  B  as  well  as  the  set  (T,pQ). 
This  aliowa  the  computation  of  extremals  for  various  g-levels, 
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The  accessory  equations  are  system  (2.26),  modified  for  the  constrained 
subarc,  and  the  discontinuities  (2.  27)  are  included  for  corner  points.  An 
additional  column  solution  is  obtained,  for  partials  with  respect  to  B,  by 
integrating  the  nonhomogeneous  accessory  equations 


(2.  77) 


where  the  no.-.homogeneous  partials  in  (2.  77)  are  computed  through  the 
appearance  of  R  in  the  control  function  (2.  73)  from  the  last  of  (2.  7l).  The 
forcing  terms  are  zero,  of  course,  for  the  unconstrained  subarc,  and  the  ini¬ 
tial  condition  for  (2.  77)  ia  the  zero  vector. 


Let  r^(t)  and  C^(t),  i=  1, ....  4;  j=  1, ....  5,  be  the  elements  of  the  partial 
derivative  solutions  obtained  from  (2.26)  and  (2.  77).  Then  the  Newton-Raphson 
equations  for  the  expanded  set  (2.26)  and  (2.  76)  may  be  written 
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2.  Experiences  with  Various  Optimization  Schemes 


The  objective  of  the  studies  was  to  obtain  a  10-g  optimal  trajectory,  and  to 
evaluate  various  optimization  schemes  during  the  process.  The  stating 
trajectory  was  the  unconstrained  optimal  trajectory  with  terminal  range 
X3  =  1000  miles,  shown  in  Figures  2  S  through  2-13,  with  peak  acceleration 
of  about  20.  7  g's.  Existence  of  the  10-g  trajectory  was  assumed  from  the 
Uj  =  16  degrees  acceleration  ploi  of  Figure  2-2.  The  function  f  of  Appendix 
B  was  taken  as 

f  =  Wj2  jvm-xj2  +  w22|?(t)  -x2/r)2  +  w32  (c(t>-x3J  2 

(2.  79) 

+  w42|Pj<t>|!  *  *52h2+«62  (B-B0)2 

The  problem,  as  stated,  turned  out  to  be  extremely  difficult,  and  a  10-g 
optimal  trajectory  was  never  found.  Thus,  the  optimization  schemes  were 
compared  under  worst  conditions. 

The  modified  Newton-Raphson  scheme,  used  first,  refused  to  converge  after 
several  successful  iterations.  Marquardt’s  method  (Appendix  B)  produced 
good  initial  changes  in  f,  and  gradually  deteriorated  in  performance  as  B 
approached  Bq.  The  difficulty  appeared  to  be  that  second-order  terms, 
neglected  in  the  development  of  the  scheme  became  larger  as  the  minimum 
was  approached.  When  the  weights  in  (2.  79)  were  changed,  the  scheme  again 
worked  well  at  firs',  then  deteriorated  rapidly.  It  always,  however,  produced 
some  improvement  in  f. 

The  parameter  B  was  reduced  to  about  11,  5  g’s  during  the  course  ol  the 
experiments,  whereas  the  other  components  in  (2.  79)  gradually  became  larger. 
The  constrained  subarc  exhibited  a  variety  of  behavior.  It  sometimes  con¬ 
tained  a  corner,  and  other  times  did  not,  and  in  some  trajectories  it  was 
absent. 


The  Fletcher-Powell  method  was  never  successfully  used.  The  H-raatrix 
(see  Appendix  B),  modified  after  the  first  step,  was  either  indefinite  (it 
should  he  positive-definite),  or  it  had  a  very  small  eigenvalue.  Thid  was  on 
the  order  of  computer  roundoff  error  and  caused  the  method  to  break  down. 

!Vo  means  of  correcting  the  situation  was  found,  although  a  recent  note 
(Reference  10)  may  contain  the  answer. 

The  use  of  the  identity  matrix  as  the  initial  H-matrix  led  to  the  optimal  gradient 
method,  in  which  the  gradient  direction  is  used,  and  the  magnitude  of  change 
is  chosen  to  minimize  f.  Initial  convergence  was  good  but  rapidly  went  to 
zero  due  to  the  "saw-tooth"  effect.  This  effect,  noted  by  many  users  of 
gradient  techniques,  makes  level  surfaces  of  f  appear  as  extremely  elongated 
ellipses,  and  gradient  changes  in  general,  therefore,  produce  little  or  no 
change  in  the  functional  to  be  minimized. 


The  use  of 


H 


o 


b<t>' 
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(2.  30) 


as  the  initial  H-matrix  (see  Appendix  B)  results  in  the  Newton-Raphson 
direction  for  the  first  step  of  the  Fletcher-Powell  method.  This  led  to  the 
optimal  Newton-Raphson  method,  which  was  tried  on  a  trajectory  which  satis¬ 
fied  B=  Bo  but  whose  terminal  conditions  were  very  far  from  those  desired. 
The  method  worked  quite  well,  and  almost  produced  the  10g-optimal  trajec¬ 
tory.  The  singular  point  u  =  ±rr,  however,  caused  computational  difficulties, 
so  the  approach  was  abandoned. 


The  predictor  scheme  of  Appendix  A  was  also  tried  on  this  problem.  In  terms 
of  the  notation  of  Equation  (2.  78),  the  derivatives  used  in  the  predictor 
equation  (2.  29)  are 
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<25(t> 
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0 

fj(0) 

f2(0) 

f3(0) 

f4<0) 

0 

(2.  81) 


A  17,  9-g  optimal  trajectory  was  obtained  using  the  optimal  Newton-Raphson 
method  (Equations  (2.  28)  with  C  determined  to  minimize  (2.  79)  with  w^,  =  0). 
The  initial  conditions  for  the  iterations  were  taken  as  those  for  the  Uj  =  55 
degrees  optimal  trajectory  of  Figures  2-1  through  2-8.  A  second  optimal 
trajectory  for  B  =  16.  9  g's  was  similarly  obtained,  using  the  17.  9-g  optima) 
values  as  starting  conditions.  A  cubic  was  fit  through  the  initial  conditions 
and  their  derivatives  ^obtained  from  (2.81)]  to  obtain  estimated  initial  conditions 
for  the  B  =  17.  7-,  17.  5-,  17.  3-,  and  17.  1  -g  optimal  paths.  The  new  paths 
were  then  obtained,  and  the  predictor  scheme  was  used  ^with  h=  0.  2  g  in  (2,29)j 
to  compute  optimals  down  to  B  =  14,  9  g’s.  Al!  of  these  paths  had  a  corner 
point,  which  migrated  toward  the  initial  point  of  the  constrained  subarc  as  B 
was  lowered.  The  predictor  scheme  failed  to  obtain  the  14.7-g  optimal  be¬ 
cause  the  corner  point  was  too  close  to  the  initial  point  of  the  constrained 
subarc.  This  caused  computational  problems  because  of  the  singular  point  u=  0. 
On  the  other  hand,  the  terminal  range  of  the  1G.  7-g  optimal  trajectory  was 
extended  to  1460  miles  using  the  predictor  scheme,  and  the  corner  point 
migrated  toward  the  endpoint  of  the  constrained  subarc.  Once  again,  the 
singular  point  u=  0  prevented  further  range  extension.  Thus,  two  parameters 
were  found  which  controlled  the  location  of  the  corner  point,  and  a  10-g  opti¬ 
mal  trajectory  could  have  been  obtained  by  manipulating  them  properly.  This 
was  not  done,  however,  since  a  lC-g  optimal  trajectory  with  a  continuous 
control  function  was  sought.  Instead,  the  constraint  equation  was  modified 
so  as  to  isolate  the  singular  point  u=  0  and  to  thus  allow  the  corner  point  to 
move  across  the  constrained-unconstrained  subarc  junction  points. 
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3.  Constraint  Equation  Modification 

The  singular  point  u  =  ±  n  has  a  clear  physical  significance:  in  the  next 
instant  of  time  the  pilot's  acceleration  will  exceed  B  g's,  since  no  control 
function  exists  which  will  alleviate  the  situation.  Computations  are  accordingly 
stopped  at  this  point.  (Normally,  0  <  |  u  |  <tt,  on  a  constrained  subarc,  and 
I  u  |  moves,  in  time,  toward  zero  or  tt  .  The  subarc  either  terminates  normally, 
or  runs  into  computational  difficulties  as  a  singular  point  is  approached.)  The 
point  u=  0  should  define  the  endpoint  of  the  constrained  subarc,  since  no  control 
exists,  for  the  next  instant  of  time,  which  will  keep  ap  =  B  g's.  It  should  be 
possible  to  show,  mathematically,  that  this  is  so.  However,  the  mechaniza¬ 
tion  of  any  such  solution  would  undoubtedly  add  complexity  to  an  already  com¬ 
plicated  computer  program.  It  was  decided,  instead,  to  modify  the  constraint 
equation  so  as  to  isolate  the  singular  point  u=  0.  The  modified  equation  reads 

a  -  „sPv2 
P  2  mgo 

With  the  constant  2K  -  0.0001,  the  added  term  normally  contributes  very  little 
to  the  original  equation;  hence,  there  is  a  very  small  difference  between  the 
two  expressions.  The  error  maintains  the  "real"  acceleration  less  than  B  g’s. 

As  u  approaches  zero,  however,  the  adder  becomes  large,  and  effectively 
isolates  the  point  u=  0.  The  equation  to  be  solved  for  the  control  function, 
analogous  to  Equation  (2.  70),  is 

a  cos2  u  +  2b  cos  u  +  c  +  r~~ ■ —  -  0  (2.  83) 

1-cos  u 

with  a,  b  and  c  defined  by  (2.  71).  Although  this  may  be  solved  as  a  cubic, 
it  was  found  more  convenient  to  use  Newton's  method.  Also,  since  cos  u  may 
be  close  to  unity,  the  transformation 

z  =  cos  u  -  1  -  •  <2. 84) 

was  introduced  for  somewhat  better  accuracy. 


7G 


The  transformed  version  of  (2.  83)  reads 

2  +  2(a+b)  +  (a+jb+c)  _  2K  _  0  (2.85) 

a  a  az 

The  constrained  subarc  equations  in  the  computer  program  were  changed 
for  the  modification.  Lack  of  time  prevented  more  than  one  or  two  debugging 
runs  to  be  made,  so  the  corner  point  was  never  removed  from  the  constrained 
subarc. 


F.  EXTENSION  TO  THE  BOUNDED  STATE  COORDINATE  PROBLEM 

1.  Statement  of  the  Problem 

The  theory  for  the  problem  of  minimizing 
T 

J  =  J  fD(x..  u)d  T 
0 


(2.86) 


subject  to  differential  equations 

x  -  f(x.  u),  x(0)  =  xQ,  (2.87) 


inequality  constraints 

G(x,  u)  2  0.  (2.88) 

in  which  u  must  appear  explicitly,  and  terminal  conditions  in  either  the  form 


'  : 

*i<T>  ’  ::i  = 

0,  i=l .  rsr. 

(2.89) 

or 

__ 

x  <T)  -  * 

J  J 

0,  j=  1,  —  ,  r- 1  s  n 

-  ; 

-  - 

(2,90) 

o' 

ii 
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i 

b 
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77 


is  treated  in  Reference  1.  The  objective  here  is  to  add  inequality  constraints 
of  the  form 

G(x)  a  0  (2.91) 

and  to  develop  a  method  for  numerically  solving  the  resulting  problem. 


2.  Necessary  Conditions 


A  set  of  necessary  conditions  for  this  problem  has  been  known  for  some  time 
(see  References  11-14).  They  are  stated  here  with  one  control  function  and  one 
inequality  constraint  equation  (2.91)  assumed  for  simplicity.  Generalization 
to  more  control  functions  and  more  inequality  constraints  is  readily  accomplished. 
It  is  also  assumed  that  the  time  derivative  of  (2.91), 


(2.  92) 


contains  the  control  function  explicitly.  The  case  where  higher  derivatives 
are  required  to  involve  the  control  function  is  treated  in  References  11  and  14. 


The  Hamiltonian  for  this  problem  is 

H  =  fQ  +  p'f  ,  (2.  93) 

where  f  is  the  integrand  of  (2.  80),  p  is  an  n  dimensional  multiplier  vector, 

(  '  )  represents  transpose,  and  f  is  the  right  hand  side  of  the  vector  differential 
equation  (2.  37).  Now  a  constrained  subarc  is  one  over  which  ineouality  (2.  91) 
is  an  equation.  A  necessary  and  sufficient  condition  for  G  tc  be  zero  over  such 
a  subarc  is  that  G  be  identically  zero  over  that  arc  (see  Reference  13).  This 
condition  is  included  in  a  new  Hamiltonian,  defined  by 

Hj  =  H  +  pG,  (2.94) 
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where  the  new  multiplier  >1  is  identically  zero  over  unconstrained  subarcs 
(G(x)  >  0)  but  may  be  different  from  zero  over  constrained  subarcs.  Equa¬ 
tions  (2.  93)  and  (2.  24)  clearly  have  the  same  value,  and  the  last  term  of  (2.  94) 
may  be  though  of  as  a  constraint  on  the  Hamiltonian. 

The  usual  necessary  conditions  nosv  ho)d.  The  Euler-Lagrange  equations 
are  derived  from  the  Hamiltonian  (2.  94),  and  the  minimum  principle  requires 
that 

H(x,  u,  p,  p)  £  H(x,  U,  p,  (J.),  (2.95) 

where  u  is  the  extremal  value  and  U  is  any  admissible  control  (satisfying 
GaO  over  constrained  subarcs).  The  Clebsch  condition  must  also  be  satisfied, 
and  usually  is  if  (2.25)  holds.  In  fact,  it  can  be  shown  that  the  Clebsch  condition 
gives  the  necessary  condition 

n  £  0  .  (2.96) 

Boundary  conditions  at  the  initial  and  terminal  points  of  the  trajectory  are 
the  same  as  those  gi\  en  in  Reference  1.  Additional  conditions  may,  however, 
be  required  at  junction  points  between  constrained  and  unconstrained  subarcs. 

If  there  are  only  two  subarcs,  the  condition  G=  0  may  be  treated  as  either  an 
initial  condition  or  terminal  condition,  depending  on  the  ordering  of  constrained 
and  unconstrained  subarcs,  and  the  multipliers  will  be  continuous  over  the 
path.  If  the  ordering  is  constrained-unconstrained-constrained  for  a  three- 
subarc-path,  the  multipliers  will  be  continuous  for  the  same  reason.  All 
other  cases  with  three  or  more  subarcs  will  produce  discontinuous  multipliers. 

It  is  well  known  (see  Reference  12)  that  the  discontinuities  take  place  at  one 
end  of  the  constrained  subarc  and  that  the  multipliers  are  continuous  at  the 
other  end.  It  f  oes  not  matter  which  end  has  the  discontinuities  (for  proof  see 
appendix  G),  so  the  initial  point  t^  is  chosen  here,  as  in  Figure  2-33. 
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Figure  2-33.  State  Coordinate  Geometry  in  the  Vicinity 
of  the  Constrained  Subarc  Initial  Point 


The  necessary  conditions  at  t^  (the  analog  of  the  Weierstrass-Erdman 
corner  conditions)  then  read 


P+<V 


P  (tj)  -  v 


3G<t}) 

dx 


(2.97) 


H 


+ 


H 


(2.98)* 


G'ttj)  =  0.  (2.99) 

where  superscripts  +  and  -  indicate  limits  from  the  right  and  left, 
respectively. 

Notice  that  when  (2.  97)  is  substituted  into  the  left-hand  side  of  Equation  (2.  98), 

•  -f. 

the  coefficient  of  v  becomes  G  (  Equation  (2.  92)]  ,  which  is  zero  by  defini¬ 
tion.  Equation  (2.98)  is  thus  independent  of  v,  and  contains  only  p  values  of 
the  multipliers.  This  equation  can  usually  be  reduced  to  an  equivalent  nec¬ 
essary  condition. 


If  G  contains  t  explicitly,  then  the  Hamiltonian  is  discontinuous  by  .he 
3G 

amount  v  -gp  (t^) , 
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For  example,  it  is  used  in  the  bounded  brachistochrone  problem  (subsection 
F4)  to  show  that  the  optimal  control  is  continuous  in  time  at  the  point  t^, 
which  implies  that 

G  '  =  0,  <2.  100) 

This  condition,  rather  than  (2.98),  is  used  there. 


3.  Basis  for  Computational  Scheme 

In  the  optimization  problem  treated  in  previous  subsections,  the  extremal 
solutions  were  functions  of  the  independent  variable  (call  it  t)  and  the  multi¬ 
plier  vector  pQ.  The  constant  v  in  equation  (2.  97)  cannot  be  determined 
from  the  necessary  ccnditions,  and  hence,  becomes  an  additional  parameter 
for  the  solutions.  The  extremal  solutions  thus  have  the  functional  forms 


x  =  x(t,  po,  v) 
P  =  P(t.  PQ,  v). 


(2.  101) 


Each  time  the  multipliers  are  discontinuous  another  constant  v  is  introduced. 
Unless  one  has  an  a  priori  knowledge  of  the  number  of  constrained  subarcs,  the 
problem  could  have  a  variable  number  of  variables .  This  gives  no  theoretical 
difficulty,  but  the  practical  bookkeeping  problems  in  a  digital  computer  pro¬ 
gram  could  become  unmanageable.  In  what  follows,  then,  it  is  assumed  that 
the  optimal  path  consists  of  three  subarcs,  ordered  unconstrained -cons train ed- 
uncons  trained. 


The  necessary  conditions  at  the  terminal  point  give  (n+1)  equations  in  the 

(n+2)  variables  (T,  pQ,  v).  The  necessary  conditions  (2.98)  (or  equivalent) 

and  (2.99)  determine  the  point  tj  and  give  an  additional  equation  in  (T,  pc,  v).  - 


The  problem  is  thus  one  of  determining  the  solution  of  (n+2)  aquations  in  (n+2) 
unknowns  and  the  Newton-Raphson  method  may  oe  used  to  find  the  solution. 

Partiul  derivatives  of  the  solutions  with  respect  to  pQ  are  obtained  as  before, 
by  integrating  the  accessory  differential  equations..  An  additional  column  in 
the  solution  matrix  is  reserved  for  partials  with  resDect  to  v.  Initial  condi  ¬ 
tions  for  this  solution  are  obtained  by  differentiating  Equation  (2.97),  and 
noting  t"at  x(tj)  is  independent  of  v.  The  method  of  oolving  the  problem  is 
illustrated  in  the  following  example. 


4.  The  Bounded  Brachistochrone  Problem 


The  bounded  brachistochrone  problem  was  chosen  to  test  the  theory.  It  ig 
simple  enough  to  have  an  analytical  solution,  yet  nonlinear  in  nature.  The 
equations  of  motion,  as  given  in  Reference  1),  are 


x  *  v  cos  y 

(2.  102) 

y  ■  v  sin  y 

(2.  103) 

v  B  -g  sin  y. 

(2.  104) 

problem  is  to  minimize  the  time  it  tukes  to  go  from  a  given  initial  point 
-  x.j  while  satisfying  the  path  Constraint 

G  ■  y-ox-biO. 

(2.  105) 

state  /ector  has  the  components  (x,  y,  v)  and  y  is 
eturmined.  Figure  2-  34  shows  the  path  constraint, 

the  control  function  to 

terminal  condition  and 

a  possible  path. 
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Figure  2-34.  Bounded  Brachistochrone  Problem  Geometry 
The  time  derivative  of  (2.  105)  is 

O  »  v(uin  y  -  u  coa  y)  (2.  lOfi) 

no  the  Hamiltonian  is  written  as 

Hj  *  1  +  PjV  cos  y  +  p2v  sin  y- p3g  sin  y +(iv(ain  y- a  cos  y)  ,  (2,107) 

The  Fulcr-Lugrange  equations  become 

Pj  •  0  <2.  tOB) 

P2  "  (2.  100) 

p3  *  - [pj  cos  y  +  p^  sin  y  ]  (2, 1 10) 

0  ■  -PjV  siny  4  p^v  coe  >  -  p^g  cos  y  +  Kvicos  v  4  a  sin  y)  <  (2,111) 
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Over  unconstrained  subarcs  the  multiplier  p  is  zero,  so  (2.  Ill)  may  be  solved 
for  y(v,  pr  p2<  p3)  as 


Pov  * 

tan  y  ■  — * - 


Plv 


or 


sin  y  = 


-(P2V  *  P3g) 


VlPlV|2  +  (P2V  '  p3gf 


cos  v  ■ 


-Plv 


Vivi2  +  i 


p2v  '  p3g 


(2.112) 


(2.113) 


(2. 114) 


The  minus  signs  In  (2.  113)  and  (2.  114}  are  chosen  to  minimize  the  Hamiltonian 
(2.  107). 

f  ** 

On  constrained  subarcs,  solution  of  G  ■  0  Equation  (2.  106) j  gives 

ton  y  *  a  ,  (2. 115) 

since  by  assumption  v  /  0,  and  Equation (2.  1 1 1) gives 


P 


P,va  -  (p2v  -  p3g) 
v(l+aZ) 


(2.116) 


Now  let  <t>  be  the  angle  defined  by  (2.  113)  and  (2.  114),  This  is  a  convenient 
definition,  since  at  the  boundary  point  t2  (Bee  Figure  2-34)  everything  is 
continuous,  so  «  becomes  y  at  that  point.  The  constrained  subarc  Hamil¬ 
tonian  may  then  be  written 

H  -  1  -  V(PlVf  +  lp2v'p3g)'!  co8<y-*>  •  (2.117) 


64 


The  minimizing  value  of  y  makes  the  cosine  positive,  so  it  must  satisfy 

*  -  <  v  <  ffl  •  4  ■  <2-  118> 

thus  resolving  the  ambiguity  of  Equation  (2.  115). 

The  multiplier  p  of  Equation  (2.  116)  is  negative  over  the  constrained  subarc. 

To  see  this,  consider  the  complete  statement  of  the  minimum  principle  for 
this  problem.  It  reads 

U(y>  *  H(  T )  (2.  119) 

for  all  T  satisfying  G(F)  2  0.  Written  out, 

PjV  cos  y  +  (p2v-p3g)  sin  ys  PjV  cos  T+  (p2v-p3g)  sin  T  ,  (2.  120) 

or  from  (2.  1  17), 

COS  <y  -  ®)  5:  COS  (r  -  (*)  .  (2.  121) 

Now  r=  <t  contadicts  (2.  121 )  if  y  4  <h ,  so  <2  must  satisfy  G(  <t> )  <0  except  at 
the  endpoint  t2  where  y  -  <£  .  Thus, 

G(  <Z> )  =  v(sintf-acos<f)<0,tj<t<t2.  (2. 122) 

Substitution  of  (2.  113)  and  (2.  114)  then  gives 

- a  gassa—  [PjV  a  -  (p2V-p3g)|  <  0,  tj<  t  <  (2,  123) 

Vlplv|2  +|P2V'P3g)2 

The  bracketed  term  of  (2.  123)  is  the  only  term  which  can  be  negative,  and  is 
the  numerator  of  (2.  1 16).  Sinc°  the  denominator  of  (2.  116)  is  positive,  p  must 
be  negative  over  the  second  subarc.  An  interesting  interpretation  of  the  angle  0 
is  that  it  is  the  extremal  value  of  y  if  the  constraint  is  not  present. 
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The  relevant  boundary  conditions  for  this  problem  are  the  terminal  and 
intermediate  conditions 


x(T,  -  xf  =  0  (2.  124) 

y(tj)-ax(tj)-b«0  (2, 125) 

and  the  transversality  conditions 

p_(T)  =  0  (2.  126) 

P3(T)  =  0  (2.127) 

PjUj)  =  p'ttj)  +  va  (2.128) 

p2^1>  "  p2^tl>  *  v  ^2‘  12£^ 

Pjttj)  *  PgUj)  (2.130) 

H+(t  j )  =  H"(tj)  (2.131) 

H  *  0  (0  st  sT).  (2.  132) 


The  condition  (2.  131)  implies  continuity  of  the  control  function  at  t=  t^,  which 
in  turn  implies  G  (t,)  =  0,  or 

v(t^)  sin  y  (t^)  -  a  cos  y  (t^)  =  0.  (2.  133) 


Continuity  of  y  follows  from  substitution  of  (2.  128)  through  (2.  130)  into 
(2.  131),  giving 


PjV(cos  y  -  cos  y  )  +  (p^v-p^gMsin  y  -  sin  y  )  +  vv^acoey 


+  + 1 
'  -  sin  y  J 


}-  °- 

(2. 134) 


where  the  parameter  tj  has  been  suppressed  and  the  bracketed  term  vanishes. 


Now  let  4/be  the  angle  defined  by  (2.  113)  and  (2.  114)  using  p  values.  (2.  134) 
then  becomes 

cos  (  y+  -  < p)  =  cos  ( y  -  =  1. 

since  y~  -  so  it  follows  that 

y+  =  4/  -  y  ,  (0  £  y  <  2  n  )  (2.135) 

on  the  optimal  path. 

From  a  numerical  standpoint  the  problem  falls  into  two  categories,  determined 
by  the  number  of  subarcs  contained  in  the  path. 

Case  1:  G  >  0  Over  the  Entire  Path  (One  Subarc)  --  Although  the  optimal  path 
is  to  consist  of  three  subarcs,  it  is  poscible  that  some  intermediate  extremal 
will  not  reach  the  boundary.  The  reduced  differential  equations  of  the  extremals 
are  (2.  1 02> - (2 .  104)  and  (2.  1 08) - (2.  1 10)  with  control  y  determined  from  (2.  1 13) 
and  (2.  114).  The  problem  then  becomes  one  of  finding  (T,  pQ)  which  satisfy 


x(T, pQ)  -  Xj  =  0  (2. 13G) 

P2  ■  0  (2.  137) 

o 

P3(T.  p0)  -  0  (2. 138) 

H(p0)  -  0.  (2.  139) 


The  modified  Newton-Raphson  equations  are 
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0 

1 

x(T) 

nn(T) 

P2 

0 

0 

1 

0 

0 

p3(T> 

= 

P3(T) 

C31(T) 

C  32^  * 

<33(T> 

H<Po> 

0 

f  1(0) 

f2<0) 

f3(C) 

_  _ 

_ 

_ 

<2.  140) 


where  the  r|'s  are  the  elements  of  the  first  row  of  r — 1 — and  the  Q  's  are  the 

5d(T)  °Po 

second  and  third  row  components  of  — -  .  Since  the  Hamiltonian  is  zero 

®  Pq 

over  the  entire  path  it  can  be  evaluated  at  t  =  0  where  all  p's  are  initial 
values,  so  the  last  row  contains  the  partials  of  the  Hamiltonian  with  respect 
to  the  pQ's  at  that  point.  Equation  (2.  136)  is  used  as  the  stopping-condition  for 
the  integrations  which  accounts  for  the  zero  in  the  left-hand  vector  of  (2.  140). 
The  accessory  differential  equations  are 


c. 


afl 

9fl 

3y 

3y 

Jy 

3  y 

a 

+  5T 

Tv 

^3  + 
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C  2  + 

i  n 
a 

1 

**2 

3f2 

o> 

3f,  | 

3y 

3y 

3y 

+  W 

3  v 

n3  +  • 

3y  j 

.8pi 

Ci  + 

3p2 

C  2  + 

dp3 

*3 

T7 

7h  ' 

3  + 

*f3 

TT 

*[£ 
i _ 

•  V 

^2 

^3 

^3 

(2.  141) 
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and  the  initial  conditions  are 


0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

1 

0 

0 

0 

0 

1 

0 

0 

0 

0 

1 

0  | 

(2.  142) 


The  ldSt  column  of  (2.  142)  will,  of  course,  give  a  zero  vector  solution 
of  <2.  141).  It  is  included  for  the  three-subarc-case  in  which  case  it  will 
contain  the  partials  of  x  and  p  with  respect  to  v. 

Case  II:  Two  or  Three  Subarcs  --  The  first  and  third  (if  present)  subarcs 
satisfy  the  differential  equations  for  Case  I.  Over  the  second  subarc  the 
reduced  differential  equations  of  the  extremals  are  (2.  102)-  (2.  104)  and 
(2.  1 00 > - (2 ,  110)  with  y  determined  from  (2.  115)  and  (2.  118).  Since  y  is 
a  constant,  the  accessory  equations  simplify  to 

Tlj  =  COS  y  Tl3 

Ti  =  sin  7  ti 


(2. 143) 


[  COS  7  Cj  +  Sin  yC  2 


The  initial  conditions  for  (2. 108)-(2.  110)  are  given  by  (2.  1 28 ) - (2 .  130),  with 
constant  v  to  be  determined,  whereas  initial  conditions  for  (2,  143)  are 


+ 

ij 


+  (xi  '  xx+)  ^  J=1*  2’  3 


i4 


*2] 


0,  i  =  1,  2,  3 

C,” 


^2 


>  j  »  1.  2,  3 


at 


C3j  s  +  (P3  *  P3>  *5p~ 


1  ) 


o 

J 


'14 


C2+4 


'34 


=  0. 


(2.  144) 


(2.  145) 


(2.  146) 


(2.  147) 


In  these  equations  j  designates  the  column  number  in  the  matrices 

4~-  *  .  and  i  corresponds  to  the  numbering  in  (2.  143). 

vPo  °Po 

Intermediate  extremals  computed  during  iterations  will  not  in  general  satisfy 
the  necessary  conditions  at  the  point  tj.  In  particular,  the  control  function 
may  be  discontinuous,  which  implies  G"  i  0.  In  this  case  the  partials  of 
in  (2.  144)  and  (2.  146)  are  computed  from  Gttj)  ■  0,  i.  e. , 


2 


1 uSill2!. 

-ax'  +  y* 


=  1,2,3  fG*  i  0). 


(2.148) 


If  G~  happens  to  be  zero,  the  coefficients  of  ^ —  disappear  in  (2. 144) 

op°j 

and  (2.  146).  These  partial  derivative  solutions  are  then  continuous  at 

t  =  t  ,  and  ^tl  need  not  be  computed.  All  variables  are  continuous  at 

3p 

*o 

J 

the  junction  point  t ^  between  the  second  and  third  subarcs  (if  there  happens 
to  be  a  third  subarc). 

Two  additional  necessary  conditions  are  required  to  set  up  the  modified 
Newton-Raphson  equations.  These  are 

G(tj)  =  y(t j ,  pQ)  -  ax(tj,  pQ)  -  b  =0 
G'ftj)  =  y"(tr  pQ)  -  ax-  (tr  pQ)  =  0  . 


It  G'(tj)  i  0,  then  (2.  149)  determines  tj(p0),  and  (2.  150)  is  the  additional 
necessary  condition  to  be  satisfied.  The  modified  Newton -Raphson  equations 
are  then 


(2.  149) 

(2. 150) 


-C 


where 
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(2. 151) 


(2.152) 
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(2.  153) 


-  o  -  2 

G(tj)  =  ag (sin  y  -  cos  V  ) , 

3t . 

and  ls  defined  by  Equation  (2.  148).  If  (2.  150)  is  satisfied, 

P°i 

then  it  determines  tj  ,  and  (2.  149)  becomes  the  additional  equation  which 
must  be  satisfied.  In  this  case  the  last  equation  of  (2.  151)  is  replaced  by 

dG(t.)  dG  SG 

-c°v  ■  —  dv*rd‘V*rdP3o 


where 


&G(t.) 

- -  =  r,„.  -  an.  .  ,  i  =  1,  2,  3 

dp  Vi  il 

*o. 
i 


The  First  Constrained  Path  --  It  is  possible  to  choose  the  initial  values 

(p  ,  p  ,  p.,  )  such  that  the  path  strikes  the  boundary  (2. 149)  before  the 
Jo  <Jo  ao 

stopping-condition  (2.  136)  is  satisfied,  'mere  will  then  be  a  second  subarc, 
and  the  problem  is  that  of  determining  a  value  of  vto  go  with  it.  Since  the 
optimal  path  is  to  consist  of  three  subarcs,  the  scheme  chosen  is  based  upon 
leaving  the  constrained  subarc  as  soon  as  possible. 


According  to  Appendix  G,  the  multiplier  discontinuities  may  be  computed 
at  either  end  of  a  constrained  subarc.  Assuming  they  are  determined  at  the 
point  t  =  t^.  Equation  (2.  110)  takes  the  form 


P3  =  -  (  p3  cos  v  +  p2  sin  y  ]  . 

Since  p7  -  p  p'  =  p9  ,  and  y+  is  found  from  (2.  115)  and  (2.  118), 
P3(f)  is  will  determined.  The  conditions  p  (t  0)  =  0,  (5  >  0  and  v  con¬ 

tinuous  must  be  satisfied  if  t^  is  the  endpoint  of  the  second  subarc.  The 
first  condition  is  solved  for  v  giving 


(p.;v  -  Pgg)  -  p-jva 
v(l+a2) 


* 

p  (t^)  i*  0  is  satisfied,  so  is  well  determined. 


(2.  154) 


(2.  155) 


<2.  156) 
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whereas  the  second  requires  that 


G  (t2)  - 


g'! 


2.  2 

-a  )  cos 


V  >  0 


(2.  157) 


Inequality  (2.  157)  is  always  satisfied  for  the  value  of  a  considered.  The 
determining  factor  is  the  continuity  of  y  which  requires  {from  the  minimum 
principle) 


sgn^Pj  +  va)  =  agn  (cos  y+) 

4. 

sgnj  (p2  ■  v)v  ■  P3S]  =  ~sen  ^sin  y  )• 


(2.  158) 


Satisfaction  of  (2.  158)  then  determines  the  endpoint  of  the  second  subarc. 

If  the  third  subarc  comes  back  to  the  boundary,  the  second  subarc  is  ex¬ 
tended  to  this  point,  and  testing  for  the  end  of  the  second  subarc  is  resumed. 

On  the  other  hand,  if  the  second  subarc  extends  to  the  endpoint,  then  ti14(T)  = 
C34(T)  =  0,  and  C24(T)  =  1.  Solution  of  (2.  151)  (or  its  equivalent  if 
c~  =  0)  gives  iterative  corrections  which  tend  to  satisfy  the  last  three  nec¬ 
essary  conditions.  The  computed  dv  is  ignored,  since  the  necessary  condi¬ 
tion  for  p^(T)  on  a  two  subarc  optimal  is  p2(T)  =  v,  rather  than  zero.  The 
optimal  path  is  to  have  three  subarcs,  so  the  corrections  are  used  iteratively 
with  the  first  constrained  path  logic.  A  three-subarc  path  is  eventually  obtained. 

Numerical  Results  --  The  constants  for  the  pi  oblem  solved  are  the  same  as 
those  used  by  Dreyfus  in  Reference  11.  They  are  xq  =  0,  y^  =  6,  v  =  1, 
g  =  32.  172,  a  =  -0.  5,  b  -  5  and  x^  =  6.  A  set  of  initial  multipliers  were 
found  which  produced  a  three-subarc  path  having  a  terminal  time  of  1.  25429 
seconds.  The  optimal  path  of  Figure  2-35  resulted  after  33  iterations,  and  a 
comparison  of  this  path  with  the  exact  solution  is  given  in  able  2-4. 
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Tabic  2-4.  Comparison  of  Computed  and  Exact  Optima!  Paths 


Quantity 

Computed 

Value 

Exact 

Value 

7  (degrees) 

-85. 2577 

-85. 2578 

0. 385117 

0.  385141 

X(tj) 

1 . G0687 

1. 60710 

y(tj) 

4.  19657 

4. 19645 

G(tj) 

1  x  10‘5 

0.  25  x  10'9 

G'(tj) 

-0, 7  x  10'3 

0 

l2 

0.533530 

0.  533526 

\(t2> 

3.  18470 

3.  18465 

y(t2) 

3.  40766 

3.  40767 

T 

0.  742246 

0.  742245 

s(T) 

6.  00001 

6 

v(T) 

2. 75571 

2.  7557] 

H(0) 

0.  1  x  10'10 

0 

Hdj) 

0.  6  x  10'15 

0 

p2m 

0. C  X  io"13 

0 

MT) 

0.7  x  10'11 

0 

Tile  differences  between  exact  and  computed  values  resulted  from  rather 
loose  interpolation  procedures  in  the  computer  program.  Although  these 
could  have  been  corrected  to  obtain  more  accurate  results,  it  was  not  felt 
necessary  since  the  purpose  of  solving  the  problem  was  to  prove  the  method. 


SECTION  III 

THE  CONTROL  FUNCTION 


A.  JUSTIFICAT'ON  OF  THE  FORM  OF  THE  CONTROL  FUNCTION 

The  form  of  the  nonlinear  optimal  feedback  control  function  is 

u  *  u(x,  t),  (3.  1) 

where  x  is  an  n-dimensional  state  vector  and  t  is  the  independent  variable. 

This  form  is  justified  theoretically,  and  a  simple  example  is  given.  The 
general  procedure  for  obtaining  the  control  is  also  outlined. 


1,  Theory 

Consider  a  set  of  optimal  trajectories,  starting  at  various  initial  conditions, 
which  all  satisfy  the  same  set  of  terminal  transversality  conditions.  The 
control  functions  for  their  trajectories  are  clearly  those  required  for  the  non¬ 
linear  optimal  feedback  control  scheme.  The  facl  that  the  control  fox  the 
family  of  trajectories  is  of  the  form  (3.  i)  follows  from  the  field  theory  of  the 
calculus  of  variations.  For  the  control  problem,  a  field  is  a  region  of  (t,  x)  - 
space  with  which  there  is  associated  a  set  of  slope  functions,  control  functions 
and  multipliers 

x  •  x  (t,  x) 

u  »  u  (t,  x)  (3(  2) 

p  "  p  (t,  x), 

having  continuous  first  partial  derivatives,  which  satisfy  the  differential 
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equations  x  =  f  (x,  v),  and  which  make  the  Hilbert  integral 


J|h(x, 


u  (t,  x),  p  (t,  x))  dt  -  p'(t,x)  dx 


independent  of  path.  A  field  is  constructed  through  the  use  of  a  theorem 
which  states  that  an  n-parameter  family  of  trajectories  which  smoothly 
covers  a  region  of  (t,  x)  -  space  and  cuts  a  surface  transversally  exactly 
once  defines  a  field.  (The  proof  fellows  from  Bliss'  corollary  83.  1  and 


the  fact  that  I*  ■  0  on  a  transversal  surface.)  Smoothness  is  established  by 

the  non-singularity  of  a  matrix  representing  >  as  in  Section  IIC,  over 

ox0 

each  path.  The  transversal  surface  is  excluded  from  the  field  since,  by 


construction,  the  matrix 


Mil) 


is  singular  at  the  endpoint  of  each  trajectory. 


T  -  K  «  0 
x2(T)  -  X2  “  0  . 


(3.4) 


"(3.5) 


(3.6) 
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The  Hamiltonian  is 


2 

H  =  ~2  +  Pix2  +  P2U' 

so  the  Euler-Lagrange  equations  are 
Pi  =  0 
p2  =  "Pi 

and  the  control  satisfies 


The  solutions  after  (3.  9)  is  inserted  into  (3.  5),  are 

t2  3 

X1  =  X10+X20t-P20T  +Plot 

t2 

x2  3  x20  “  P201  +  P10  2 

Pi  =  PfO 

p2  =  P20  ~  PlO1, 

where  tne  (0)  subscripts  indicate  conditions  at  t  ■  0.  The  transversality 
conditions  (Reference  1)  require 


Pl(T>  »  0 

P2(T)  =  v2 
H  =  v3. 


(3.7) 


(3.8) 


(3.9) 


(3.  10) 


(3.11) 


which  is  non-singular  everywhere,  except  at  t=K.  The  trajectories  thus 
smoothly  cover  the  solution  space.  The  first  two  of  (3.  13)  may  be  solved  to 
obtain 


100 


(3.  16) 
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so 'that  (3.  19)  becomes 


x  -  x(t,xo) 

P  =  P<t,xo). 


(3.21) 


Note  that  the  transversality  conditions  and  terminal  equations  must 
constitute  an  n-dimensional  non-singular  (or  normal)  system  to  obtain  the 
solution  (3.  20).  Non-singularity  of  the  matrix  (3.  15)  allows  inversion  of 
the  first  of  (3.  21)  to  obtain 

x  =  x  (t,  x) ,  (3.  22) 

o  o 

which  is  the  general  form  of  (3.  1G).  Substitution  of  (3.  22)  into  the  last 
of  (3.21)  gives 

p  =  p(t,  x)  (3.  23) 

as  the  general  form  of  (3.  17).  The  control  (3.  18)  then  assumes  the 
desired  form  (3.  1).  Substitution  of  the  control  into  the  original  differential 
equations  then  gives  the  last  of  the  field  equations  (3.  2). 


In  more  complex  problems,  the  process  described  above  can  be  done  only 
in  a  small  region  about  a  trajectory,  as  in  Section  II  C. 


B.  METHODS  FOR  DETERMINING  THE  CONTROL 

As  shown  in  the  example,  the  control  may  sometimes  be  obtained  by  direct 
solution  of  the  optimization  problem.  This  is  generally  impossible;  however, 
one  might  include  in  the  category  of  "direct  solution,  the  faster  than  real 
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time  solution  of  the  optimization  problem  in  an  airborne  computer.  This 
approach  is  unrealistic  at  the  present  time,  due  to  large  computational 
requirements. 

An  alternate  possibility  is  that  of  solving  of  an  approximate  optimization 
problem.  The  model  differential  equations  are  usually  as  simple  as  the 
process  allows;  however,  most  optimal  trajectories  are  fairly  insensitive 
to  path  deviations,  so  far  as  the  optimal  criterion  is  concerned.  With  feed¬ 
back,  the  approximate  problem  solution  would  approach  an  optimal  solution 
as  the  re-entry  progressed  so  that  terminal  errors  would  be  nulled  as  well. 

Some  time  was  spent  seeking  suitable  simplifications  to  the  re-entry  problem, 
and  some  useful  relationships  were  found.  However,  the  problem  was  not 
realistically  simplified  to  the  point  that  analytical  solutions  could  be  obtained. 
Another  possibility  considered  was  that  of  using  a  Rayleigh-Ritz-Galerkin 
procedure  to  compute  near-optimal  trajectories  in  the  airborne  computer. 
This  process  also  ran  into  large  computer  requirements,  and  hence,  was 
not  pursued  further. 

A  final  alternative  is  that  of  using  surface  fitting  procedures  to  approximate 
the  optimal  control  function.  This  has  the  advantage  of  moderate  on-board 
computer  requirements  and  may  well  handle  large  deviations  from  "nominal" 
re-entry  conditions.  The  method  also  has  growth  potential,  in  that  self 
adaptive  features  may  be  built  in  by  increasing  the  dimensionality  of  the  fit. 
However,  several  facts  should  be  borne  in  mind  in  such  an  approach.  There 
is,  in  general,  no  way  of  determining  the  best  form  of  the  approximating 
function.  Polynomial  approximations  are  probably  best  from  the  standpoint 
of  evaluation  in  the  airborne  computer,  so  this  form  is  assumed  for  the  re¬ 
mainder  of  the  report. 


1C3 


The  region  over  which  an  approximation  is  to  be  valid  is  usually  determined 
by  the  errors  which  can  be  tolerated.  It  is  generally  not  true  that  higher- 
ordered  fits  produce  more  accuracy  after  a  certain  point,  for  computer 
truncation  and  roundoff  errors  destroy  the  benefits  of  the  additional  terms. 

It  is  better  to  segment  the  region  into  subregions  and  to  use  lower-ordered 
fits  to  obtain  more  accuracy.  The  segmentation  is  readily  accomplished 
and  mechanized  for  one -dimensional  approximations,  but  may  lead  to  serious 
mechanization  problems  in  more  dimensions. 

The  larger  the  dimensionality  of  the  approximation,  the  harder  the  problem  of 
obtaining  the  fit.  It  is  best,  if  possible,  to  split  the  over- all  approximation 
problem  into  subproblems  of  reduced  dimension. 

The  polynomial  approximation  approach  was  chosen  as  the  approach  for 
further  development  of  the  non-linear  optimal  feedback  control  scheme. 


C.  DEVELOPMENT  OF  THE  CONTROL  FUNCTION 


The  polynomial  approximations  considered  are  of  the  form 


w 


.  i  j  k  1 

V-t/l  *2  *3  2  ' 


i*  j.k.  I 


i.j.k  ,1  20, 
i+j+k  srn, 


(3.  24) 


where  w  is  the  approximated  control  function  u,  or  one  of  the  multipliers 
Pj,  P2 .  z  is  range,  and  the  variables  y^y£  anc*  ^3  are  rested  toV,  y  h 
through  the  equations  .  _ 


*1 

=  v-vn 

*2 

c 

1 

i'¬ 

ll 

(3.  25) 

^3 

=  h-hn  . 
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Nominal  state  variables  V  ,  y  ,  hn  arc  taken  as  those  for  trajectory  1 
in  Table  2.  2,  and  the  reason  for  the  transfor  mation  is  made  clear  below. 
The  multiplier  approximations  are  used  in  evaluating  the  expression 


u  =  tan 


-  1 


~CLOP2 
2t  DLP1V 


(3.26) 


It  is  possible  to  split  the  over-all  approximation  (3.  24)  into  two  parts 
because  of  the  method  in  which  the  data  was  collected  (Section  II  C).  Thus, 
consider  the  two  approximations 


Vz)  ■ 


b. ..  ,2 

ijk  l 


,1*0 


Va,,._(z)  y2*  y3k  .  i. j.kaO 


i.  j.k 


i+j+k  sm 


The  method  of  solution  is  to  determine  sets  of  coefficients  a.^(  zr), 
zr  =  15r  ,  r  =  0,  1,  ....  99, 

for  the  three-dimensional  approximation  (3.  28).  and  to  use  these  as 
data  for  the  several  one-dimensional  fits  (3.27). 


Consider,  then,  the  problem  of  determining  the  coefficients  a^j.  for 
a  given  value  of  nnge,  It  is  convenient  to  transform  the  problem  to 
the  form 


w 


\  ilk 

£>  >  2  3 

ijk 


i.  j.k  20 
i+j+k  sm 


(3.27) 


(3.28) 


(3.29) 


(3.  30) 
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si,  n  =  1,2,3,  through  the  relationships 


in  which  I  x 
'  n 


xi  =  °i  >i  +  n 

x2  =  a2  y2  +  >2.  (3.  31) 

x3  =  a3  y3  +  y3* 

I 

The  transformation  (3.  3l>  gives  a  better  conditioned  system  of  equations 
to  be  solved  for  the  coefficients,  and,  moreover;  makes  it  easy  to  find 
negligible  coefficients,  if  there  should  happen  to  be  any.  The  coefficients 
a.  ..  are  obtained  by  substituting  (3.  31)  into  (3.  30)  and  rearranging  the 
resuits. 


When  derivatives  are  included  as  data,  tKe  normal  procedure  is  to 
differentiate  the  approximating  polynomial  to  obtain  equations.  Thus,  for 
example. 


3w  dw  3V  dy  1  3w 

dXj  "  dV  dyjdxj  dV 


Y  Cijk  ixll  x2’’x3k* 

ijk 


with  similar  expressions  for  ^  .  The  four  eouations 


dx- 


(3.  32) 


- 

2 

2 

2 

1 

X1 

X2 

X3 

X1 

xlx2 

xlx3 

x2 

X2X3 

X3 

w 

0 

1 

0 

0 

2X1 

X2 

x3 

0 

0 

0  •  •  • 

c  = 

dw 

^xj 

0 

0 

1 

0 

0 

X1 

0 

2x2 

x3 

0  •  •  • 

Sx2 

0 

0 

C 

1 

0 

0 

X1 

0 

x2 

2x3  ’  '  * 

j 

#0(^> 
_ 1 

(3.  33) 


with  vector 


C100’ 


c010’  C00T  c200’  C1 10’ 


c101’  c020’  C0ll’  C002’ 


IOC 


ar  e  expansions  of  (3.  30)  and  its  partial  derivative  equations,  in  which 
terms  through  second-order  in  x  are  shown.  The  matrix  of  (3.  33)  may  be 
evaluated  numerically  for  each  of  the  27  trajectories.  These  are 
collected  into  a  matrix  system  of  equations 

At  =  f 

in  which  A  is  a  (108  x  q)  matrix,  c  is  a  q- dimensional  vector  of  co¬ 
efficients,  and  f  is  the  right-hand  side  of  (3.  33)  for  the  27  trajectories. 
The  dimension  q  is  related  to  the  order  of  fit  m  through 

q  =  1  +  ( 1 1  +6m  -t-m  ^) , 

and  the  range  on  m  is  0<m  «C. 


1 .  Lagrange  Fits 

The  approximating  lunction  is  required  to  agree  with  actual  function 
values  and  derivatives  q  times.  Then  q  of  Equations  (3.  34)  are  selected 
and  the  right-hand  sides  contain  the  relevant  function  values  and  partial 
derivatives  determined  during  the  data  generating  process.  The  question 
of  which  equations  to  US'*  is  resolved  or.  the  basic  of  the  least  singular 
system  of  equations,  as  determined  by  the  process  of  Appendix  F.  Notice 
that  the  matrix  ot  evaluated  basis  functions  is  the  same,  regardless  of 
which  function  (u,  Pj,  P la  fit. 

An  exploratory  program  was  written  to  obtain  Lagrange  fits  for  (u,P.,P.,) 
for  all  values  of  m  considered.  As  originally  constructed,  the  matrix  of 
bu^ls  functions  was  inverted,  and  three  matrix  multiplies  produced  the 
Coefficients  for  the  fits,  Much  fit  was  evaluated  at  the  27  data  points,  and 
errors  were  computed.  Output  included  actual  and  computed  function  values 
errors,  and  the  menu  and  deviation  for  each  fit,  computed  from 


(3.34) 


(3.  35) 
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E 


(3.  36) 


27 


.y  ii 

27 


o 


(3.  37) 


in  which  is  the  error  in  the  function  for  the  i**1  trajectory.  Additionally, 
the  arc  tangent  relationship  (3.26)  was  evaluated,  and  errors  and  atatis.'ics 
for  this  were  included  as  output.  Adequate  precautions  were  taken  to  ensure 
th^t  the  inverse  matrix  was  well  conditioned.  If  the  product  of  the  matrix 
of  basis  functions  and  its  inverse  wore  not  close  to  the  identity  matrix,  a 
well-known  Newton-Raphson  correction  was  applied.  *  Those  function  errors 
which  were  supposed  to  be  zero  provided  an  additional  accuracy  check  (at 
least  cnc  function  equation  must  be  included  in  the  set  of  equations  solved, 
to  obtain  the  coefficient  c00q). 

The  first  attempts  at  fitting  used  ( V,y,  h)  in  place  of  (Xj,  x2>  Xg)  in  Equation 
(3.  30).  The  matrix  of  basis  functions  turned  out  to  be  very  badly  conditioned, 
so  the  linear  transformations 


V 

“1 

5  »1  V 

X2 

=  «2  y 

+  02 

- 

(3.  38) 

X3 

=  e-3  h 

+  ^3 

with  |  xi  |  s  l,  1=1, 2, 3,  were  introduced.  The  cr's  and  0's,  of  course,  are 
computed  from  maximum  and  minimum  values  of  (V,  y,  h)  at  the  given  value 
of  rang<n  The  new  matrix  of  basis  functions  was  well  conditioned,  and  fits 
for  all  values  of  m  were  obtained. 

*Sce  Reference  9,  report  No.  3,  pg.  216.  The  NASA  project  reports  contain 
a  wealth  of  information  on  surface  fitting  techniques,  as  applied  to  a 
problem  similar  to  the  one  considered  here. 


The  exploratory  program  was  run  for  several  range  points  and  the 
following  conclusions  were  reached: 

•  Fits  for  m  =  5  in  Equation  (3.  30)  were  beBt  for  most  of  the 
range  points  explored.  Maximum  control  errors  on  the  order 
of  10°  were  found,  although  the  means  and  deviations  were 
much  less. 

•  The  control  function  fit  was  superior  (in  all  cases)  to  the 
fit  obtained  from  Equation  (3.  25). 

•  The  equations  used  in  the  fitting  process  changed  from  one 
range  point  to  the  next. 

•  The  equations  used  in  the  fitting  process  changed  when  the 
accuracy  of  the  constants  in  Equation  (3.38)  was  increased 
from  5  to  8  significant  digits  (at  a  given  range  point),  in¬ 
dicating  great  sensitivity. 

•  The  constants  in  Equation  (3.  38)  were  dependent  on  range. 


An  attempt  was  made  to  fit  the  constants  of  Equation  (3.38)  se  functions 
of  range  so  that  the  approximations  (3,30)  could  be  used  directly.  The 
ChebyHcheff  fitting  program  described  below  was  used  in  the  attempt. 


In  view  of  the  behavior  of  the  variables  (see  Figures  2-16  -  2  28)  it  is 
evident  that  the  constants  are  only  piecewise  smooth  as  functions  of 
range.  Hence,  the  fits  obtained  were  not  good.  Moreover,  when  the 
fitted  constants  were  used  in  the  transformation,  the  inverse  basis 
function  matrix  again  became  ill  conditioned.  Hence,  it  became  neces¬ 
sary  to  perform  a  "back"  transformation  by  inserting  (3.38)  into  (3.  30) 
and  rearranging  the  results.  This  gave  polynomials  in  terms  of  (V,  y,  h). 

The  coefficients  a.  .  so  obtained  were  fairly  large,  indicating  that  subtraction 
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must  take  place  in  the  evaluation  of  the  polynomial.  In  a  seven  digit  word- 
length  computer,  this  could  lead  to  severe  truncation  and  roundoff  error. 


Moreover,  the  back  transformation  could  not  be  performed  at  all  for  fits 
like  (3. 30)  in  which  3  i  i,  j,  k  a  0,  (The  computer  used  had  a  12-digit  word- 
length,  and  the  fits  obtained  were  somewhat  worse  than  those  for  ma  5 
above.)  When  the  transformation  (3.31)  was  used,  the  back  transformation 
was  readily  accomplished,  and  the  coefficients  a...  so  obtained  were  wall 
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behaved.  Thus  <he  conclusion: 


•  The  approximation  (3.  28)  with  variables  (3. 1.5)  should  be  used  for 
evaluation  in  a  seven  digit  wordlength  airborne  computer. 

Tne  transformation  (3.  31)  is  obtained  from  (3.  38)  by  adding  and  subtracting 
.he  appropriate  amounts.  Thus, 


X1  =  =  o,  (V-V)  r(0lv  +1'31)  =  a1y1+y1 


(3.39) 


with  similar  expressions  for  X2  and  x^. 

It  was  decided  to  explore  least  squares  approximations  as  the  next  step  in 
the  development  of  the  control  function. 


2.  Least  Squares  Fits 


One  may  readily  verify  that  the  system  of  equations  to  be  solved  is  (see 
Reference  2) 


(A  A)  c  =■  A 


(3.40) 


no 


where  the  notation  is  that  of  Equation  (3.  34),  and  vector  f  contains  function 
values  and  partial  derivatives  obtained  during  the  data  generating  process. 

Unit  weights  are  used  since  the  derivative  equations  were  of  the  same  order 
of  magnitude  as  the  functional  equations.  An  exploratory  program  similar 
to  that  for  the  Lagrange  fits  was  constructed,  except  that  a  subroutine  to 
solve  simultaneous  equations  was  used  in  place  of  matrix  inversion.  The 
accuracy  of  the  fits  so  obtained  was  very  good,  as  measured  by  the  closeness 
to  zero  of  the  mean,  Equation  (3.  36).  One  may  readily  verify  that  the  first 
normal  equation  for  the  system  (3.  34)  requires  (3.  36)  to  be  zero.  Then  the 
deviation.  Equation  (3.  37),  may  be  interpreted  as  the  functional  portion  of 
the  error  minimized  by  the  least  squares  process. 

The  exploratory  program  was  run  with  the  following  conclusions: 

•  The  least  squares  fits  were  better  than  the  corresponding 
Lagrange  fits. 

•  Fits  for  m  =  5  were  again  best  for  most  of  the  range  pointB. 

•  The  control  function  fit  was  superior  (in  all  cases)  to  the  fit 
obtained  from  Equation  (3.  26). 

•  The  approximation  (3.  28)  with  variables  (3.  25)  Bhould  be  used 
for  evaluation  in  a  seven-digit  wordlength  airborne  computer. 

It  was  decided  to  use  the  least  squares  approximations  for  the  control  function. 
Accordingly,  fifth-order  least  squares  fits  were  obtained  for  the  range  values 
(3.29),  excluding  r  ■  95  -  99.  The  order  of  the  best  fits  for  these  range 
points  arc  summarized  in  Table  3-1. 


Ill 


Table  3-1.  Order  of  Beet  Fite  Near  Terminal  Point 


Range  Point 

95 

96 

97 

98 

99 


Order  of  Fit 

(m) 

4 

4 

1 

1 

2 


The  control  functions  of  Figure  2-28  are  shown  in  Figures  3rl  through  3-27, 
together  with  the  corresponding  errors  in  the  least  squares  approximations. 

The  maximum  errors  occur  in  the  vicinity  of  maximum  change  in  the  control 
function  which  corresponds  to  the  bottom  of  the  first  dip  into  the  atmosphere. 

The  approximations  are  quite  good  before  this  region,  and  very  accurate  from 
about  z  =  800  to  z  =  1300  miles,  which  covers  the  complete  skip  after  the  first 
dive.  The  approximations  become  somewhat  worse  as  the  endpoint  is  approached. 


*?.  Chebyscheff  Fits 

The  final  phase  of  the  control  function  approximation  is  that  of  performing  the 
56  one -dimensional  approximations  (3.27).  This  was  not  done,  since  time  did 
not  permit  both  this  approximation  and  the  simulation  studies  of  Section  IV  to 
be  accomplished.  It  is  evident,  from  Figure  3-20,  that  the  one -dimensional 
fits  must  be  segmented  into  several  parts  in  order  to  obtain  any  degree  of 
accuracy  in  the  approximations  so  the  final  phase  is  not  a  trivial  problem. 

A  Chebyscheff  criterion  fitting  computer  program  was  constructed  for  this 
final  phase.  The  method  is  due  to  Professor  F.  Koehler  of  the  University  of 
Minnesota,  who  was  a  consultant  during  the  contract  period.  Assume  that 
equally -spaced  data  is  given  over  the  interval  ^0,  lj  ,  and  that  the  order  of  the 
fit  to  be  obtained  is  n.  A  Lagrange  polynomial  of  order  n  is  constructed  for 
the  (n+1)  tabular  points  closest  to  the  Chebyscheff  zeros,  denoted  by 
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Figure  3-2.  Optimal  Control  u  and  Error  of  Fit  for  Trajectory  2 


Figure  3-3.  Optimal  Control  u  and  Error  of  Fit  for  Trajectory  3 


Figure  3-5.  Optimal  Control  u  and  Error  of  Fit  for  Trajectory  5 


Optimal  Control  u  and  Error  of  Fit  for  Trajectory  6 


Figu 


Control  u  and  Error  of  Fit  for  Trajectory  19 


Control  u  and  Error  of  Fit  for  Trajectory  20 


Figure  3-21.  Optimal  Control  u  and  Error  of  Fit  for  Trajectory  21 


Control  u  and  Error  of  Fit  for  Trajectory  22 


Figure  3-24.  Optimal  Control  u  and  Error  of  Fit  for  Trajectory  24 


Control  u  and  Error  of  Fit  for  Trajecto 


coefficient  aQ00  (radians) 


Max  I  f(x)  -  Pn<x)  !  =  Min .  (3. 42) 

[o  s  x  si] 

In  general,  this  is  not  true,  the  method  proceeds  iteratively  to  the  solution. 

Let  e(x)  be  the  error 


e(x)  =  f(x)  -  Pn<x)  (3. 43) 

which  vanishes  at  n+1  points  by  construction.  If  the  data  is  smooth  enough, 
e(x)  will  vanish  exactly  n+1  times,  and  there  will  be  n+2  maximum  errors 
(two  at  0  and  1,  and  n  relative  maximum  and  minimumsj.  The  objective,  on 


each  iteration,  is  to  equalize  the  n+2  absolute  errors.  Let  z o>  Zj, . . . ,  *n+1 
denote  the  points  at  which  the  maximum  errors  occur,  and  let  q(x)  be  a 
correction  polynomial  chosen  to  equalize  the  errors.  Then 

-e(zo)+  h^)  =  e  (lijl-ntej)  =  ...  3  (-l)n  [c<zn+1?  " 


Let.  h  be  the  common  value  in  (3.44)  so  that 


n(zQ)  =  h  +  c(zQ)  - 

^(•Zj)  a  *h  +  e  Uj)  . 

(3.45) 

Tl(z  )  •>  (-l)“h  +  c(z  ) 

n  n 

1<W'  +  E<zn+1> . 
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hquotions  (3.45)  are  n+2  equations  in  h  and  the  n+1  coefficients  of  the 
nl  ’  order  polynomial  n(x).  The  approximating  polynomial  for  the  next 
iteration  is  taken  as 

Pn(2)(x)  PnU)(x)  +  h(x).  (3.4 

Iteration  is  required,  since  the  location  of  the  maximum  errors,  in  general, 
shifts  at  each  stage. 

The  method  was  found  to  work  very  well  in  practice.  It  rapidly  converged 
to  stationary  points,  z.,  i=  0,  1,  .  .  .  ,  n+l,  and  the  selection  of  the  same  set 
of  points  z.  on  two  successive  iterations  was  found  to  be  a  good  stopping 
condition. 


SECTION  IV 


MECHANIZATION  AND  SIMULATION  OF  THE 
NONLINEAR  OPTIMAL  FEEDBACK  CONTROL  SCHEME 


SIMULATION  DESCRIPTION 


The  mechanization  of  the  schene  for  the  re-entry  problem  is  illustrated  in 
Figure  4-1.  Detailed  descriptions  of  each  block  in  the  diagram  follow. 


INITIAL 
CONDITIONS  x 


Figure  4-1.  Nonlinear  Optimal  Feedback  Control  Scheme 
Mechanization  for  the  Re-entry  Application 

1.  Re-entry  Vehicle  and  Sensors 


In  view  of  the  short  time  of  flight  (15  minutes  or  less),  accelerometers  are  the 
only  sensors  assumed  to  be  aboard  the  vehicle.  They  measure  the  acceleration 
components  (smoothed,  if  necessary)  in  the  lift  and  drag  directions. 


a,  *  -^2Cir  sinu 
1  2m  LC 

a2  =  ‘TB“(CDO  +  CDL  C0S  Uj  • 
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The  constants  ,  C^q,  C£>q  and  ^DL  ma^  be  c^'an8ec*  from  their  nominal 
values  to  evaluate  the  over-all  performance  of  the  scheme,  and  the  density 
profile  may  assume  one  of  several  forms  for  the  same  purpose. 


2.  Navigator 


The  navigator  integrates  the  differential  equations  of  motion  and  supplies  the 
state  vector  and  its  time  derivative  at  discrete  instants  of  time.  In  terms  of 
accelerometer  measurements,  the  equations  cf  motion  are 


2 

dv  gQR  sin  y 

=  a  o - o — 

dt  (R  +  hr 

2 

dy  v  cos  y  gQR  cos  Y 

dt  v';  (R+h)  o(R+h)2 

dh 

=  v  sin  y 

dz  R  v  cos  y 

at  =  Cg«+HT  • 


(4.2) 


The  state  vector  components  are  velocity  v,  flight  path  angle  Y ,  altitude  h, 
and  range  z  (measured  in  miles),  and  the  state  vector  derivative  is  the  right- 
hand  side  of  (4.  2).  The  total  heat  is  also  found,  for  scheme  evaluation  purposes, 
by  integrating  equation  (2.  2)  along  with  system  (4.  2), 


3.  Predictor  and  Control  Generator 


The  predictor  equation  is  the  Adams-Moulton  equation 


n+1 


=  x 


A 


55x  * 
n 


59  x 


n-i 


+  37x 


n-2 


9x 


n-3 


(4.3) 


used  to  introduce  lead  into  the  system.  Lack  of  time  prevented  experiments  with 
other  predictor  equations,  and  A  =  1/2  and  2  seconds  were  the  only  time  lucre 

ments  considered. 
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Equation  (4.3)  requires  three  back  derivatives  as  well  as  the  present  point, 
so  it  is  not  applicable  for  the  first  6  seconds  of  flight.  The  first  three  back 
derivatives  are  generated  from  (4.  2)  with  a  constant  control  in  Equations 
(4.  1).  The  actual  value  of  this  control  is  unimportant  since  the  veMcle  is 
out  of  the  sensible  atmosphere  during  the  first  6  seconds.  However,  it  is 
taken  as  the  control  generator  equation  value  at  the  initial  point  (which 
assumes  zq  -  0)  for  continuity  purposes.  (Figure  2-28  shows  that  th 
control  is  almost  constant  during  the  initial  phase  of  the  flight.  > 


In  Section  111,  100  control  function  fits  of  the  form 


u 


i.J.k  *0 
i+J+k  £  m. 


(4.4) 


in  which  V^,  and  h^  are  "nominal"  trajectory  values,,  were  found.  The 
coefficients  a. ..  were  noted  to  be  functions  of  range;  however,  they  were  not 
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reduced  to  polynomial  form  due  to  lack  of  time. 


Tne  control  generator  interpolates  linearly  between  fitted  range  points  (spaced 
15  miles  apart)  to  determine  the  coefficients  a^  at  the  predicted  range  point. 
Also,  the  nominal  trajectory  values  are  determined  by  the  same  process.  The 
polynomial  (4.  4)  is  then  evaluated  with  the  predicted  values  of  V,  Y  and  h. 


Of  course,  the  control  is  constant  for  the  first  6  seconds  of  flight.  It  is  also 
evaluated  differently  for  the  last  30  miles.  This  is  necessitated  by  the  singu¬ 
lar  endpoint  and  the  general  compi  ession  of  the  state  vector  differences  near 
the  endpoint  (see  Figures  2-16  through  2-28).  Hence,  a  linear  interpolation 
for  the  control  is  performed  with  range  as  the  variable.  The  two  points  used 
are  the  control  at  zn+j  -  1470,  and  -nat  z  =  1500  miles,  z  ■  1500  miles  is 
the  stopping  condition  for  simulated  trajectories. 
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4 .  Hold  Circuit 


The  hold  circuit  supplies  the  constant  control  during  the  first  6  seconds,  and  the 
linear  function  of  range  during  the  last  30  miles  of  flight.  In  the  intermediate 
region,  it  is  taken  as  a  cubic  in  time  using  the  points  u  ,  un_2  to  evaluate 

the  coefficients.  The  transformation 


S 


2  + 


(t  -  v 

A 


(4.  5) 


allows  the  cubic  coefficients  to  be  computed  from 


(4.  6) 


B.  SIMULATION  RESULTS 


Table  4-1  contains  results  of  simulation  runs  for  the  27  optimal  trajectories. 

AJ  is  the  difference  between  simulation  and  optimal  trajectory  total  heat.,  and  is 
seen  to  be  a  positive  quantity;  A  is  the  time  increment  used  in  the jAdams- 
Moulton  predictor  equation  (4.  3),  and  Av  and  Ah  are  the  terminal  velocity  and 
altitude  errors  from  the  optimal  trajectory  values  of  1650  ft/sec  and  75,  000  feet 
respectively. 


Table  4-1.  Simulation  Results  for  the  27  Optimal 
Trajectories  of  Table  2-2. 


Trajectory 

Number 

Total  Heat 
Optimal  Traj. 

Total  Heat 
Simulation 

AJ 

(BTU/ft'H 

Ah 

(ft) 

1 

28, 872 

28, 898 

26 

2 

226 

-  8, 655 

2 

31, 888 

31, 913 

25 

2 

59 

-  8,  565 

3 

26, 252 

26, 270 

18 

2 

234 

-  8,  537 

4 

29, 278 

29, 295 

17 

2 

218 

-  8,  982 

5 

28. 666 

28, 721 

55 

2 

298 

-  5,804 

6 

29, 032 

29, 051 

19 

2 

227 

-  8,  898 

7 

28, 757 

28, 810 

53 

2 

Bl 

- 10,  566 

8 

32, 280 

32,  295 

15 

2 

226 

-  8,884 

9 

32,  116 

32,  148 

32 

1/2 

281 

-  6, 254 

10 

26,657 

26, 687 

30 

2 

211 

-  9, 018 

11* 

25, 792 

25, 844 

52 

2 

-  49 

-25,  640 

12 

32,  135 

32, 344 

209 

1/2 

220 

-  7, 649 

13 

31, 844 

31, 900 

56 

2 

246 

-  9,34..' 

14 

26, 362 

26, 382 

20 

2 

225 

-  8,  952 

15 

26,  128 

26, 173 

45 

2 

244 

-  9,  520 

16 

29, 963 

Fail  at  z  = 

1056  miles 

17 

29, 017 

Fail  at  z  = 

1431  miles 

18 

28, 510 

28, 562 

52 

2 

243 

-  9,  603 

19 

28, 817 

Fail  at  z  = 

472  miles 

20 

•>  AAA 

O  *-»  ,  i  *  * 

Fail  at  z  - 

1077  miles 

21 

31,374 

31,397 

23 

2 

249 

-  8.447 

22 

31,  1 10 

31,  173 

63 

2 

250 

-  9, 168 

23 

3 1 , 375 

31, 449 

74 

1/2 

263 

-  5,  570 

24 

27, 326 

27,348 

22 

2 

220 

-  9, 048 

25 

26, 891 

27,043 

152 

1/2 

137 

-12,  201 

26 

26, 216 

26, 266 

50 

2 

253 

-  8, 930 

27 

26, 368 

26, 417 

49 

1/2 

-  43 

-10, 750 

*  Az  =  -6  miles. 
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Most  of  the  trajectories  were  simulated  using  A  -  2  seconds.  The  advantage 
of  making  A  as  large  as  possible  is  that  the  onboard  computer  can  be  used  for 
other  tasks  (such  as  computing  display  information)  during  the  idle  time. 

It  was  found  that  several  trajectories  failed  with  this  increment  (trajectories 
9,  12,  16,  17,  19,  20,  23,  25  and  27),  and  a  failure  was  recorded  if  the  con¬ 
trol  for  a  given  trajectory  exceeded  ±500  degrees.  This  was  an  adequate  test, 
since  if  the  control  for  a  given  trajectory  exceeded  ±180  degrees,  it  also  ex¬ 
ceeded  ±500  degrees.  Computations  were  also  stopped  if  the  altitude  fell 

below  50,000  feet.  This  happened  with  trajectory  11,  which  stopped  with  a  j 

range  error  Az  =  -6  miles,  as  shown  in  Table  4-1.  j 

j 

S 

WithA=  1/2  sec. ,  only  trajectories  16.  17,  19  and  20  failed.  This  time  i 

trajectory  1  1  stopped  with  a  range  error  A  z  =  -0.  1  mile,  as  indicated  in 
Table  4-2.  Trajectory  1  was  also  run  and  only  small  improvement  in  terminal 
conditions  resulted  (Table  4-2). 

A  simulation  program  with  no  prediction  (A=  0)  was  constructed  to  test  the 

ultimate  capabilities  of  the  control  function  polynomial  approximation.  (The  j 

polynomial  is  evaluated  at  each  integration  step,  and  the  system  is  assumed  to  | 

react  instantaneously.)  This  was  run  for  trajectories  1,  9  and  11  with  the  , 

results  shown  in  Table  4-2.  There  was  no  change  in  terminal  conditions  for  j 

trajectory  1  (compared  with  the  A  =  1/2  second  case),  and  only  slight  changes  ■ 

for  trajectory  9.  Path  11  reached  the  stopoing  condition  (z  =  1500  miles)  even 
through  the  final  altitude  was  less  than  50,000  feet.  Trajectories  16,  17,  19  I 

and  20  failed  again  when  run  with  this  program.  t 


Table  4-2.  Additional  Simulation  Results 


T  rajectory 

Total  Heat 

— 

Total  Heat 

AJ 

A 

Av 

Ah 

Number 

Optimal  Traj. 

Simulation 

(BTU/ft2) 

(sec) 

(ft /sec) 

(ft) 

1 

28, 872 

28, 872 

28, 897 

28, 897 

25 

25 

1/2 

0 

K2£9ii 

-  8, 655 

-  8, 655 

9 

32,  116 

32, 147 

31 

0 

281 

-  6,  279 

11* 

25, 792 

25,  842 

50 

1/2 

-171 

-25, 107 

11 

25, 792 

25, 842 

50 

0 

-170 

-25, 043 

28 

29, 589 

Fail  at  z  = 

1044  miles 

29 

29,011 

29, 048 

37 

2 

223 

-  8,  579 

30 

28, 496 

28,  551 

55 

2 

2  53 

-  9,026 

-rf- 

31 

28, 687 

28,  720 

33 

1/2 

263 

-  8,081 

*  Az  =  -o. : 

mile 

-  .iv--. 

• 
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Four  additional  trajectories,  numbered  28-31  in  Table  4-2,  were  run  to  test 
the  nonlinear  optimal  feedback  control  scheme  for  off-design  trajectories. 

The  initial  conditions  for  these  are  given  in  Table  4-3,  and  they  correspond  to 
the  initial  conditions  of  trajectories  20-27  with  Ax ^ q  -  0.  Trajectories  29  and 
30  ran  with  A=  2  seconds,  whereas  A  =  l/2  second  was  required  for  path  31. 
Trajectory  28  fails,  even  with  A=  0. 


Table  4-3.  Initial  Conditions  for  Trajectories  28-31 


Trajectory 

Number 

Initial  Velocity 
(ft/sec) 

Initial  Flight 
Path  Angle 
(degrees) 

Initial 

Altitude 

(ft) 

Sign  of  Change 

From  Trajectory 

**10 

A*3C 

28 

35,  000 

-5.40 

430, 000 

0 

+ 

+ 

29 

35,  000 

1 

o 

370, 00C 

0 

+ 

- 

30 

35,  000 

l 

430, 000 

0 

- 

+ 

31 

35,  000 

H 

370, 000 

0 

- 

- 

Path  28  may  be  compared  with  paths  20  and  24  to  see  the  effects  of  initial 
velocity  changes,  since  Ax2Q  and  Ax3Q  are  the  same.  Numbers  20  (+  AxJ0) 
and  28  (AxJ0  =  0)  fail,  whereas  path  24  (-AXjQ)  is  successful  with  A=  2  seconds. 
Similar  comparisons  are  shown  in  Table  4-4. 


Table  4-4,  Effects  of  Initial  Velocity  Change  on  A 


Trajectory 

Number 

Sign  of 

A 

AX10 

wsa 

21 

+ 

+ 

2 

29 

0 

+ 

► 

2 

25 

- 

+ 

- 

1/2 

22 

•  + 

- 

+ 

2 

30 

0 

- 

2 

26 

- 

- 

+ 

2 

23 

+ 

- 

1/2 

31 

0 

• 

1/2 

|  27 

- 

- 

- 

1/2 
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The  characteristics  of  the  successful  flights  are  illustrated  in  Figures  4-2 
through  4-5.  Only  those  for  trajectory  1  are  displayed,  since  most  of  the 
flights  are  fairly  similar.  The  first  dip  into  the  atmosphere  is  quite  critical, 
and  the  inaccuracies  in  the  control  function  (Figure  3-1)  drastically  affect  the 
remainder  of  the  path.  The  vehicle  emerges  into  the  skip  region  (where  the 
control  function  fit  is  very  accurate)  following  a  new  optimal  path.  This  is 
pursued  chrough  the  second  dip  to  a  point  (about  1410  miles)  where  the  con¬ 
trol  fit  again  becomes  somewhat  inaccurate.  Inaccuracies  from  this  point 
onward  do  not  affect  total  heat  very  much,  but  do  contribute  to  terminal 
velocity  and  altitude  errors.  Another  sources  of  these  errors  is  the  linear 
fit  of  the  control  in  the  last  30  miles  of  flight.  The  control  is  still  rising 
at  1470  miles  (Figure  2-28)  and  the  additional  lift  available  to  the  optimal 
trajectory  is  neglected  in  the  simulation.  This  causes  the  vehicle  to  dive 
more  steeply  and  (in  general)  to  lose  less  velocity  as  the  endpoint  is  approached. 
The  terminal  errors  should  be  reduced  by  starting  the  linear  interpolation 
closer  to  the  endpoint,  however,  no  experiments  of  this  nature  were  performed. 

It  is  understandable  that  trajectory  16  failed.  This  path  was  difficult  to  obtain 
during  the  mapping  process,  indicating  that  it  was  near  the  edge  of  the  corridor, 
and  consequently,  that  the  control  fit  is  near  the  edge  of  the  region  of  appli¬ 
cability.  Control  function  inaccuracies  could  easily  move  the  vehicle  into  a 
region  where  the  control  function  fit  is  no  longer  applicable.  A  comparison  of 
optimal  and  simulated  trajectories  shows  that  the  vehicle  is  quite  far  from  the 
optimal  at  the  point  of  failure  (Av  =  -4  50  ft/ sec,  Ah  =  7, 000  feet). 

The  large  inaccuracies  in  the  control  function  of  trajectory  17  (Figure  3-17) 
occur  in  the  vicinity  of  the  first  dip  into  the  atmosphere  (the  bottom  is  at 
z  =  600  miles).  They  cause  the  skip  to  be  quite  different  from  that  of  the 
optimal  trajectory.  The  effect,  however,  is  not  noticed  until  near  the  end  of  the 
trajectory,  since  the  accurate  control  fit  produces  a  reasonable  control  for  the 
intermediate  portion  of  the  flight.  The  vehicle  eventually  leaves  the  region  of 
applicability  for  the  control  function  polynomial. 
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gure  4-2.  Simulation  Results  for  Trajectory  i  -  Altitude 
Versus  Range 


Simulation  Results  for  Trajectory  1  -  Velocit; 
Versus  Range 


Simu 

Versi 


The  initial  conditions  for  trajectory  19  cause  the  vehicle  to  dive  deeper  into 
the  atmosphere  which  in  turn  produces  greater  lift  and  drag  forces.  The 
vehicle  follows  the  optimal  path  through  the  bottom  of  the  dip  (at  z  =  375  miles'. 
However,  the  control  inaccuracy  bulge  in  the  region  z  =  400  miles  to  z  «  4o0 
miles  (Figure  3-19)  causes  the  vehicle  to  be  captured  by  the  atmosphere, 
instead  of  following  the  skip  of  the  optimal  path.  This  shows  the  sensitivity 
of  the  paths  to  the  control  function  errors  ncur  the  bottom  of  the  first  dip. 

The  flights  of  trajectories  20  and  28  are  similar  to  that  of  path  16,  and  they 
show  that  simultaneous  positive  perturbations  in  initial  flight  path  angle  and 
altitude  are  not  tolerated  hy  the  control  function  fit. 

The  conclusions  drawn  from  the  simulation  results  are: 

•  In  its  present  form,  the  nonlinear  optimal  feedback  control  scheme 
produces  reasonable  re-entry  trajectories  with  modest  increases 
in  totai  heat  from  the  optimal  values,  for  a  suitably  restricted 
region  of  initial  conditions. 

•  The  inaccuracy  of  the  control  function  fit  over  the  first  dip  into  the 
atmosphere  drastically  affects  the  remainder  of  the  flight.  The 
accuracy  of  the  fit  should  be  made  better  here  to  eliminate  the 
failures,  and  to  consequently  decrease  the  total  heat. 

•  It  is  believed  that  the  terminal  point  errors  can  be  reduced  by 
starting  the  linear  interpolation  for  the  control  as  a  function  of 
range  closer  to  the  endpoint. 

•  The  time  interva1  A,  used  in  the  predictor  equation,  was  found  to 
be  a  function  of  the  simulated  trajectory;  A  should  be  made  a  vari¬ 
able  quantity,  dependent  upon  the  state  of  the  vehicle. 


SECTION  V 

CONCLUSIONS  AND  RECOMMENDATIONS 


A.  ACCOMPLISHMENTS 

A  re-entry  optimization  problem  is  presented  and  solved  through  the  methods 
developed  in  Reference  1.  These  produce  an  accurate  solution  to  the  optimiza¬ 
tion  problem.  A  powerful  predictor  scheme  was  developed,  which  supplements 
the  work  of  Reference  1  and  which  allows  the  optimal  solutions  to  be  changed 
as  a  function  of  a  parameter.  This  was  used  to  extend  the  range  of  the  original 
optimal  trajectory,  and  to  perform  an  "absolute  minimum"  test.  The  latter 
showed  that  the  optimal  path  is  a  minimizing  path,  at  least  over  a  large  region 
of  solution  space.  Sufficiency  tests  for  a  relative  minimum  were  also  developed. 
It  was  shown  that  the  trajectories  considered  are  relative  minimums. 

The  predictor  scheme  was  used  to  map  the  re-entry  corridor,  and  27  optimal 
trajectories  spanning  the  corridor  were  obtained.  The  control  functions  for 
these  trajectories,  and  partials  of  the  control  with  respect  to  the  state  vector, 
were  used  to  obtain  a  polynomial  approximation  for  the  control  function  over 
the  optimal  re-entry  corridor.  The  approximation  was  used  in  the  mechaniza¬ 
tion  of  a  nonlinear  optimal  feedback  control  scheme.  Simulation  results  showed 
that  modest  increases  in  total  heat  over  optimal  values  were  experienced,  and 
that  large  (although  tolerable)  terminal  errors  resulted.  It  is  believed  tha*  the 
terminal  condition  errors  can  be  greatly  reduced,  if  necessary.  A  few  of  the 
trajectories  failed,  which,  in  effect,  limits  the  region  of  applicable  initial 
re-entry  conditions  to  some  extent. 

Another  re-entry  optimization  problem,  in  which  the  sensed  acceleration  is 
constrained  to  be  less  than  or  equal  to  a  given  number  of  g's,  was  posed,  and 
several  optimization  methods  were  used  in  an  attempt  to  obtain  a  10-g  optimal 
path.  The  methods  all  failed,  although  the  predictor  scheme  emerged  as  the 
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most  powerful  ox  the  group.  A  singular  point  on  the  constrained  subarc 
prevents  convergence  to  the  solution.  A  means  of  circumventing  thiB 
problem  was  found  but  remains  to  be  tested. 

The  optimization  methods  were  also  extended  to  include  the  bounded-state 
coordinate  problem. 

B.  RECOMMENDATIONS 

The  nonlinear  optimal  feedback  control  scheme  in  its  present  form  produces 
near  optimal  re-entry  trajectories  if  the  space  of  initial  conditions  is  suitab'y 
restricted.  The  scheme  is  technically  sound  but  fails  in  some  instances  only 
because  of  control  function  errors. 

The  polynomial  approximation  is  based  upon  control  values  for  several  optimal 
trajectories,  and  the  paths  are  readily  generated  through  use  of  the  predictor 
scheme.  It  is  possible  to  generate  many  more  paths  than  were  used  in  the 
study,  and  consequently,  to  fill  out  the  optimal  re-entry  corridor  with  optimal 
values  of  the  control  function.  The  increased  number  of  data  points  should 
enable  better  control  function  approximations  to  be  made. 

The  data  processing  task  is  that  of  reducing  the  data  points  to  a  form  readily 
implemented  in  the  control  scheme.  Polynomial  approximations  are  probably 
best  for  this  purpose;  however,  there  are  many  ways  of  approaching  the  fitting 
problem.  It  may  be  that  surface  fits  of  the  partial  derivatives  are  well  behaved 
over  a  large  region  of  the  corridor,  so  that  fits  of  these,  followed  by  an  inte¬ 
gration,  may  produce  more  accurate  control  function  approximations.  On  the 
other  hand,  it  may  be  necessary  to  segment  the  corridor  into  several  pieces, 
and  to  perform  surface  fits  for  each  of  these  to  obtain  the  required  accuracy. 
Also,  it  may  be  possible  to  split,  the  over -all  fit  into  several  low-dimensional 
subfits,  so  that  each  stage  of  the  fitting  process  may  be  carefully  controlled. 
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The  data  processing  task  is  seen  to  be  in  an  incomplete  state  at  the  present 
time.  Since  many  practical  optimal  feedback  control  problems  could  be 
implemented  if  this  task  were  accomplished  in  a  practical  way,  it  is  re¬ 
commended  that: 

•  Prai  ical  methods  should  be  sought  for  performing  tne  data 
processing  task. 

The  control  function  polynomial  approximation  is  presently  based  upon 
unconstrained  optimal  trajectories.  However,  many  real  optimization 
problems  necessarily  contain  constraints  upon  the  control  and  state  vectors. 

It  is  recommended  that; 

•  Mechanization  problems  associated  with  inequality  constraints 
should  be  examined. 

In  view  of  the  rising  importance  of  adaptive  control  capabilities,  the  methods 
should  be  extended  to  this  case.  Theoretically  there  is  no  problem,  although 
practically,  the  dimensionality  of  the  surface  fit  is  increased.  Thus: 

•  Problems  associated  with  extending  the  methods  to  include 
adaptive  control  capabilities  should  be  examined. 

Finally,  most  feedback  control  schemes  which  operate  over  a  finite  time 
interval  contain  a  singular  point  at  the  endpoint.  This  is  true  for  the  problem 
studied  in  the  report.  The  singular  point  causes  difficulties  in  the  vicinity  of 
the  endpoint,  and  usually  results  in  terminal  point  errors.  Thus: 

•  The  singular  terminal  point  problem  should  be  examined 
theoretically  to  determine  if  the  endpoint  errors  can  be 
minimized. 


15h 


MATERIALS  SECTION  (TABLE  LXXV) 


No  work  of  the  nature  described  in  Section  2  of  Table  LXXV  was  performed 
under  Contract  AF33(6 1 5)- 1 858,  BPS  Number  4(6399-62405364-822501). 
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APPENDIX  A 
A  PREDICTOR  SCHEME 


Suppose  that  one  has  an  n-dimensional  system  of  equations  of  the  form 

ip  (x,  b)  =  0  (A.  1) 

in  which  x  is  an  n-dimensional  vector  and  b  is  a  single  parameter.  Assume 
further  that  the  matrix  exists  and  is  non-singular  ever  the  range  of  b 
considered,  and  that  the  vector  also  exists.  Then  by  the  implicit  function 
theorem  it  follows  that 


x  =  x  (b). 


(A.  2) 


and 


dx  _  Idjp 
db  |9x  3b 


(A.  3) 


If  m  solutions  (A.  2)  to  the  system  (A.  1)  are  known  for  equally  spaced  values 
of  the  parameter  b,  and  if  the  corresponding  derivatives  are  computed  from 

(A.  3)  (call  the  complete  system  . and  V  ,  .  ,  ,  x^1  respectively), 

the  problem  is  that  of  predicting  the  next  member  of  the  family  of  solutions, 

xm  +  r 


Open  type  integration  formulas  are  well  suited  for  this  task,  and  many  such 
formulas  are  given  in  Chapter  6  of  Reference  2.  In  particular,  a  formula 
truncated  after  third  differences  (the  Adams-Moulton  predictor  equation)  is 


m  +  1 


=  x 


m 


_h_ 

24 


(55  x  '  -  59  >: 


m 


m-T 


37  x' 
m- 


O  -  9x  '  ,), 


(A.  4) 
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where  h  is  the  spacing  between  parameter  values.  Notice  that  the  present 
point  and  derivative,  and  three  previous  derivatives  are  required  for  this 
equation.  Other  formulas  using  less  back  information  are  easily  derived 
from  the  results  given  in  Reference  2.  The  simplest  ouch  equation, 
using  only  the  present  point  and  derivative,  is  the  Newton-Raphson  equation. 


x  , ,  =  x 

m  +  1  m 


+  hx' 


m 


(A.  5) 


As  an  example  of  the  use  of  the  predictor  scheme,  consider  the  range 
extension  problem  of  Section  I1B.  The  equations  corresponding  to  (A.  1) 
are 


V(T.Pq)  -  Xx  =  0 
?(T,F0)  -  X2/R  =  0 

(A.  6) 

C(T,  P0>  -  X3  =0 
P.,(T,  P  )  -  0 

c.  O 

H  (F  )  =  0  . 

o 


The  parameter  b  is  identified  with  the 
(A  3)  become 


T  ' 

V(T) 

riii(T) 

P  ' 

10 

£<T) 

^31<T) 

P20 

= 

C(T) 

t141(T) 

P30 

p2(T) 

C2i<t) 

P40 

0 

fj(0) 

terminal  range  X 

g,  so  liquations 

n12(T) 

h23(T) 

-1 

r 

o  i 

v32<t) 

r]33^Tf 

h34<T) 

0 

r42(T) 

n43m 

h44(T) 

1 

<22<T> 

C23(T) 

C24(t> 

0 

f2(t» 

f3(° ; 

f4(0) 

0 

(A-  7) 
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It  is  seen  that  the  third  column  of  the  inverse  Newton-Raphson  matrix  contains 
the  derivatives  for  the  predictor  equation  when  the  solutions  (T,  P  )  satisfy 
system  (A.  6)  (for  a  given  value  of  the  parameter  X^). 


APPENDIX  B 

THREE  OPTIMIZATION  TECHNIQUES 


In  addition  to  the  modified  Newton-Raphson  method  and  the  predictor 
scheme,  three  other  optimization  techniques  were  considered  during  the 
period  devoted  to  optimal  trajectory  computations.  They  include  the 
optimal  Newton-Raphson  scheme,  the  Fletcher-Powell  method  (Reference  3) 
and  Marquardt's  scheme  (Reference  4), 

As  a  common  basis  for  discussion,  let  the  system  of  equations  to  be  solved 
by  the  optimization  schemes  be 


M  -  0,  (B.  1) 

where  both  \p  and  x  are  n-dimensionr1  vectors.  Multiply  this  system  by 
a  diagonal  weighting  matrix  \V  to  obtain 

p(x)  =  Wi//(x)  =  0.  (B.  2) 

W  clearly  does  not  change  the  solution,  but  does  aid  in  the  numerical 
computations.  Define  a  function 


f(x)  = 


(B.  3) 


whose  absolute  minimum  value  (zero)  corresponds  to  the  solution  of 
(B.  1)  or  (B.  2),  and  let  g(x)  be  the  gradient  vector. 


g(x) 


9d(x) 
3  x 


(B.4) 
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Both  the  optimal  Newton-Raphson  and  rietcher-Powel)  methods  require 
that  a  minimum  value  for  f  be  found  in  a  direction  s.  Let 

x  =  xq  +  is  (B.  5) 

where  xQ  is  the  present  point  and  X  is  a  parameter,  and  determine  a 
point 


x,  =  x  +  Kr[S,  0  <  K  <  1, 
1  o'4  ■ 


(B.  6) 


where 


•n 


minimum  of 


(B.  7) 


Equation  (B  7)  was  fourid  in  the  appendix  of  Reference  3,  where  W.  C. 

Oavidor.  (Reference  5)  is  given  credit  for  its  origin.  The  factor  K  in 
Equation  (B.  6)  was  found  necessary  for  the  re-entry  problem  since 
Equation  (B.  7),  in  some  instances,  produced  too  large  an  estimate  for  n,. 

If  f(xQ)  g(xo>,  f(Xj)  and  g<Xj),  are  known  numerically,  then  the  function 
f(x)  may  i  e  approximated  by  the  cubic  equation  , 

f(xQ  +  \s)  =  aQ  +  a2  X  t-  a2\2  +  a3X3.  <B.  8) 

Values  for  the  coefficients  are  found  to  be 


ao 

=  f<*o> 

al 

=  g'<xo;s 

a2 

=  ‘kT1  Z4«^o)s1 

(B.  9) 

a3 

- - 4(  2Z  +  g'(x0)  s+g'(x  )s] 

3(Kn) 

Z 

=  K|[f(xo)  ■  f(x1,]+  e'(x0)s+g/(xi)s  ■ 

* . 
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When  Equation  (B.  8)  is  differentiated  and  evaluated  for  the  minimizing 
value  V  =  a,  it  is  found  that 


-a2  +'\l a2  -  3aja3 


(B. 10) 


so  the  approximate  minimum  value  of  f  in  the  direction  s  is  given  by 


f{\Q  +  us). 


The  direction  s  in  the  optimal  Ne  .vton-Raphson  method  is  given  by 
the  Newton-Raphson  change 

f  9<Mx  Jl'1 

N  l  ?x  J  v  o 


(B.  11) 


When  K  and  r,  are  both  unity  in  (B.  6)  and  (B.  7),  x^  is  the  Newton- 
Ruphson  predicted  point,  so  the  method  becomes  straight  Newton-Raphson 
as  the  solution  is  approached. 

The  Fletcher -Powell  direction  is 


S,.  =  -H(x  )  g(x  ), 

1<  o  °  o 


where  ll(x  )  is  a  positive -definite  matrix  to  be  updated  after  each  step 
(Reference  3  suggests  the  identity  matrix  as  a  suitable  initial  choice). 
When  u. 


(B.  12) 


x^  +  Q  Sr.,, 

o  £ 


(B.  13; 


1  ..aid'll.,  i  .link,  kill 


f(x2)  and  g(x2)  have  been  evaluated,  the  updating  equations  for  the  next  step 
are 


H(x„)  =  H(x  )  +  A  +  B 
£  o 

A  -  - 

o'y 

(H(x  )y)(H(x  )y) ' 

O  -  °  U 

y'H(xo)y 
0  =  a  sf 

y  =  g(x2)  -  g(x0). 


(B.  14) 


In  Marquardt's  scheme  it  is  assumed  that  a  linearized  expansion  adequately 
describes  the  surface  behavior  in  a  suitably  small  sphere  of  radius  6o  about 
x  .  Thus,  Equation  (B.  2)  may  be  written 

d<6(x  ) 

<MxQ  +  6)  =  0  (xQ)  +  ax  9  6,  (B.  15) 

where  the  change  6  is  to  satisfy  the  constraint 

6 '6  -  6  *  =  0.  (B.  16) 


The  function  <B.  3),  using  approximation  (B.  15),  is  minimized  subject  to  the 
constraint  (B.  16)  by  setting  the  partials  of 


F  (6,\)  =  f(xQ  +6)  +  \(6‘  6  -  i>Q) 


(B.  17) 


to  zero.  Upon  performing  the  operations  it  turns  out  that 

r. 

d<t>  (x  J I  I  90<x„)  I  I 

+  xjl  a  - 


a\ 


c»x 


-g(xo> 


(B.  18) 
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In  practice,  X  '*  0  is  a  parameter.  When  X  =  0,  it  follows  from  (B.  4)  and 
(B.  11)  that  5  is  the  Newton-Raphson  change.  On  the  other  hand,  if  X  is 
large,  6  approaches  a  gradient  change.  The  policy  is  to  make  X  as  small 
as  possible  during  iterations. 
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APPENDIX  C 

THE  TRANSFORMED  OPTIMIZATION  PROBLEM 


The  system  of  differential  equations  for  the  re-entry  problem  are 


q 

v 

y 

i 


c 

yrr 

-s 

2  m 


S 

=  -2HT 


v_ 

R 


1/2  3 

P  v 


+  7.  5  N 


p  v2  C 


D 


(u)  - 


gp  siny 
(1  +  l)2 


v  cos  y 

pvC.  (u)  + - 

L  R  <1  4  5  ) 


sin  y 


l  3/2 

1  v  12.  5 

) 

10, ooO 

gQ  cos  y 
v(l  +C)2 


dC  _  v 

ar-  rrryy  cos^ 


The  transformation 

z  3  ~£~  •  cq  8  5-  280  ; 

3 


changes  the  units  of  range  from  feet  to  miles,  and  transforms  the  last 
of  (C.  2)  to 


dz  =  J_  d£_  a 
dt  c  ^  dt 


c3(l  A) 


cos 


y 


(C.  1) 


(C.  2) 


(C.  3) 


(C.4) 
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and  note  lhat 


dV 

dz 


(C .  8) 


Also,  transform  altitude  to  units  of  feet  through 


h  =  R£ 


(C.  9) 


Then  Equations  (C.  6)  become 


dci  _  C3C  (R+h)pl/2V 

dz  '  rVn“  cos  y 


dV 

dz 


7.  5  x  10 


c.^N 


R 


(B  +  hi 

cos  y 


3/2 


' c_i 

mR 


ft)  p V  g  /  »  -  2c  e  R 
cos  y  £c3S0n 


tan  y 
(R  +  h> 


.  (C.  10) 


dy  _  C35 

dz  ’  2 rnP. 


(R+h)p 

COS  >■ 


CL(u) 


C3gp  R 

V  (R-h) 


dh 
d  z 


c„ 

'“p-  (R  +  h)  tan  y 


dt 

dz 


3  (R  +  h) 

R 7Tf2  . 

V  cos  y 
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The  Hamiltonians  for  the  system  (C.  1)  and  (C.  2),  and  the  system  (C.  10) 
are,  respectively. 


H 


P'f 


and 

H2  =  So  +P'S- 

Each  Hamiltonian  represents  a  canonical  system  of  equations,  so  a 
canonical  transformation  may  be  used,  in  part,  to  transform  one  system 
into  the  other.  A  second  transformation  is  required  for  the  change  of 
independent  variable. 


The  generating  function  F2  (q,  P,  t)  of  Reference  d,  page  240,  represents 
a  canonical  transformation  between  old  and  new  state  vectors  (q,  Q)  and 
multipliers  (p,  P).  The  transformed  equations  satisfy 

3F2 

P  "  7^" 


Q  = 


dp 


K  =  H  + 


3F2 

at 


Let 


F 


2 


2 


■h 


K 


2 


v 


+  \3r?  +  \4 


C  3 


(C.  11) 


(C.  12) 


(C.  1,3) 


(CM) 
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Then  in  the  transformed  system. 
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where  K  is  the  multiplier  vector  (C.  16),  and  the  first  variation  of  the  function 


J  =  J  F  dt 
0 


(C.  18) 


is  to  be  taken.  Let  the  independent  variable  of  (C.  18)  be  changed  from  t  to 
x^  =  z.  This  may  be  done  since  range  is  an  increasing  function  of  time. 
Then  (C.  18)  becomes 


(C.  19) 


where  u  is  taken  as  the  time  derivative  of  a  fifth  state  coordinate.  The 
generalized  Lagrangiap  for  (C.  19)  is 


F 


1 


dx  . 

r  3 

- 1  \ 

i  n 

k  *  I  \ 

fi  '  "3F~I  +  K4  | 

| f  4  '  TR-] 

L  i=i 

l  "3T  1 

\  “ar/J 

dt 

dz 


(C.  20) 


In  the  following,  the  derivative  subscripts  represent  partials  with  respect 
to  the  derivatives.  One  may  readily  verify  the  relationships 


^dx.  s  -  1  =  1.  2>  3 

_ l 

dz 


,1 


u  dxj 

dz 


4 

■as 


F  ♦  £  x.i. 


i=l 


(C.  21) 

<C.  22) 

(C.  23) 
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In  (C.  23),  the  definition  for  the  Hamiltonian 

H  »  F-£  y/Fj,  (C.24) 

i  1 

has  been  used,  where  the  derivatives  y/  include  the  control  function 
as  well  as  the  derivatives  of  the  state  components  with  respect  to  the 
independent  variable.  In  the  present  case,  Fu  of  (C.  22)  is  zero  (from 
the  transversality  conditions)  so  that  (C.  23)  is  recognized  as  the  usual 
control  problem  Hamiltonian.  When  the  Hamiltonian  for  is  formed, 
according  to  (C.  24),  and  (C.  21),  (C .  22),  and  the  center  expression  of 
(C.23)  are  inserted,  one  finds 


H'  = 


-X 


4  • 


On  the  other  hand,  when  Hj  is  used  for  (C.  23)  one  finds 


I 

i  =  1 


dt 

cTF 


dx5 

F  +  X.f . 
dz  U  4  4 


(C.  25) 


(C.  26) 


The  last  expression  is  zero,  so  with  F^  taken  as  zero,  this  is  the  same 
as  the  Hamiltonian  (C.  12),  for  by  the  transformation  (C.  5)  (taken  as  the  last 
transformation). 


i  =  0,  1,  2,  3. 


(C.  27) 


APPENDIX  D 


Lot 


U 


V  = 


_&*.  +  _®2L  K 

dx  dp 

o  ro 


o  1  o 
Then  from  differentiation  of 


(D.  5) 


UU'1  =  I 


(D.  6) 


it  follows  that 


‘  - 1  - 1  •  - 1 
U  =  -U  U  U  . 


(D.  7) 


Substituting  from  (D.  1),  and  premultiplying  by  V,  one  finds 


V  U"1 


VU'1  +  VU'1-— vu 

dpdx  9p2 


Also  from  (D.  1),  by  post- multiplying  the  expression  for  V  by  U  *, 


V  U 


.-1 


r  b2H 


3x 


-  * 

J~  +  3xT5 


vu 


-1 


(D.  8) 


(D.  9) 


When  (D.  8)  and  (D.  9)  are  added,  the  left-hend  side  is  found  to  be  a  perfect 
differential,  so  that 


■dr'vu'1)  - 


32h 


a2n 


-1  3  H 


3x 


t  +  37*?  VU  *  VU  W  +  VU 


,-l  3  If 


3p' 


VU 


<D.  10) 


From  (0.5),  (D.4)  and  (D.3)  it  is  found  that 


VU-  .  -|f .  .  K. 


(D.  11) 


where  is  known  to  be  a  symmetric  matrix.  Other  fields  mav  he  generated 
cx  •*  J 

by  changing  the  initial  conditions  to  other  symmetric  matrices,  such  as  the 
mutrlx  considered  In  the  sufficiency  test,  Section  III). 


lUO 


APPENDIX  E 

RELATIONSHIPS  FOR  SUFFICIENCY  CONDITIONS 


NON-SINGULAR  NEWTON  -  RAPHSON  MATRIX  AND  NORMALITY 


are  end  conditions  satisfied  by  the  extremal  path,  and  repeated  subscripts 
in  (E.  2)  and  (E,  3)  denote  sums  .  If  the  order  of  abnormality  q  is  zero,  the 
path  is  normal, 
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where  (E.  9)  is  omitted  if  the  terminal  value  of  the  independent  vari.aole  is 
not  specified.  The  terminal  equations  (I).  2)  and  (E.3)  are  then 


St ( r)  "  'en+l+k  '  k  =  1 ,  ....  r 
Cj(T)  =0.  j  =  r+1,  - n 

e(n+l)+kXk  (  r)  =  ~en+r+2’ 


(E. 10) 


(E. 11) 


Assume  first  that  (E,  9)  is  included.  Then  the  first  r  and  the  last  (n-r)  of 
(E.  7)  at  t  =  T  =  K  give 


B  e  =  0  , 


(E. 12) 


where  B  is  the  Newton-Raphson  matrix  (as  in  Section  1IC).  If  B  is  non- 
sinimlar,  e-  0,  and  the  result  follows.  If  (E.  9)  is  omitted,  then  (E.ll) 
mu  it  be  included,  since  with  en+r+2  =  ^  an  additional  constraint 

equation.  However, 

c(n+i)+k\<T)  =  '  Ci<T)  x^T)  +p.(T)tli(T)  =  -^(0)^(0)=  6.^(0)=  0 

(E. 13) 


(E.  14) 

and  (E.  1)  when  differentiating  the  second  of  (E.  13)  with  respect  to  t  to  obtain 
zero.  The  Newton-Raphson  matrix  is  formed  with  the  total  differentials 


since  the  expression  is  constant  in  t.  This  may  be  shown  by  using  the 
canonical  equations 


lb3 


0  =  dx^IT)  =  xk  (T)  dT  +  (T)  k  =  1 . r 

0  =  ap.  =  p.(T)  dT  +  C  (T)  j  =  r  +  1, _ n 

J  J  J 


(E.  15) 


and  the  last  of  (E.  14),  as  in  Section  1LA.  Then  if  B  is  non-singular  the 
set  (e,  e  )  in  (E.5)  is  zero,  and  the  result  again  follows. 


RELATIONSHIP  BETWEEN  MAXIMUM  PRINCIPLE  AND  MINIMUM  PRINCIPLE 
ACCESSORY  SOLUTIONS 

The  difference  in  maximum  and  minimum  principle  formulations  of  the  opti¬ 
mization  problem  is  the  signs  of  the  multipliers,  including  the  unit  multiplier 
fov  the  integrand  f  of  J,  the  function  to  be  minimized.  Thus,  the  Hamiltonian 
for  the  maximum  principle  is  the  negative  of  the  Hamiltonian  for  the  minimum 
principle.  Let  the  functional  form  of  the  solutions  for  the  maximum  principle 
formulation  be 


x  =  x(t,  x  ,  z  ) 
o  o 

z  =  z(t,xo.zo) 

and  the  corresponding  forms  for  the  minimum  principle  formulation  be 
x  =  x(t,xo,po)  -  ;  ' 

p  =  p(t.x0,p0)  • 

It  follows  that 

x(t.xo.P0)  =  *<t,xo,zo) 

if  the  multipliers  for  the  two  systems  satisfy 
pU.Xo.Pq)  =  -2(t,xo,zo). 


(E.  16) 


(E,  17) 


<E. 18) 


(E.  18) 
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The  accessory  solutions  arc  the  partial  derivatives  of  (E.  18)  and  (E.  19)  with 
respect  to  xq,  po  and  z  .  Upon  differentiating,  the  correspondence  between 
systems  of  solutions  is  found  to  be 


-5; —  (t,x  ,p  ) 
ox  0  0 

o 


\ —  (t,  x  ,p  )  = 
op  o  *o 

o 


(t,  x  ,2  ) 
ox  0  0 


3x  .  3  2 

; -  (t,X  ,2  )  _ C 

*  2  O  O 

O  dp 


(t,  X  ,  p  )  =  - 

ox  '  o  Ko 


- (t,  X  ,2  ) 

ox  0  0 

o 


dx  ,,  . 

_ (t,  x  .  z  ) 


(E.  20) 


do  -fc 
z  „  V  _ o  _  &Z 


Z^VPo1  =  -TT^'W 


dz^(t'xo'V  5p  =ST"  ft'V*©* 


o  o 


In  constructing  a  field  of  extremals  for  the  sufficiency  test.  Bliss  uses  the 
maximum  principle  formulation  of  the  problem.  A  suitable  choice  for  the 
conjugate  system  of  solutions  is 

U  -  —  (t  x  ,  2  )  +  — 1^—  (t, x  ,2  )[  sxP  +  I  ] 

dxo  0  °  dzo  0  0  L  dxo  1 


V 

dx 

0  L 

0 

i  bz 

1  0 

o> 

l  0 

v  ;1r!t'vV  +  tt“-v*o>  [73T-  +  1j- 

o  O  1  o  J 

Upon  substitution  from  (E.20),  these  become,  for  the  minimum  principle, 

u  =  -^-<t'Vpo)  +-^r{t-xo’po}  [^r  '  l] 

o  ro  u 

-v  =  iir{t'x o'Pq)  (t-xo-po)[4^2--  !]  • 

O  O  u  o  J 


<E.  21) 


<E.  22) 
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REDUCTION  OF  FIELD  FROM  (n+2)  DIMENSIONS  TO  (n+1)  DIMENSIONS 


Consider  the  system  of  equations 
x  =  4r(x'z'zn+l) 


n+1 


u(x'z'zn+l> 


(E.  23) 


-z  = 


n+1 


=  0 


which  arises  from  the  maximum  principle  formulation  of  the  optimization 
problem  when  (the  integral  of  u)  is  considered  as  an  additional  state 

coordinate.  It  is  tacitly  assumed  that  xr+1  does  not  appear  in  the  right-hand 
side  of  system  (E.  23)  and  that  it  does  not  appear  in  the  boundary  conditions. 
Hence,  as  a  necessary  condition. 


n+1 


=  0. 


(E.  24) 


and  the  first  and  third  of  (E.  23)  are  the  system  of  equations  to  be  solved  for  the 
optimal  trajectory.  It  is  not  immediately  obvious  that  the  field  of  extremals 
used  in  the  sufficiency  teat  can  be  correspondingly  reduced  from  (n+2)  dimensions 
to  (n+1)  dimensions.  Tins,  however,  can  be  shown  as  follows.  The  accessory 
equations  corresponding  to  (E.  23)  are 


d2H 

0 

I  a2H 

32H  ] 

«*> 

^zdx 

!  dz2 

^zn+l 

] 

M 

d 

•3F 

Vi 

3u  1 
TT 
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du 1 
TT 
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3  z  , 
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32H 

0 

3u 

c 

'T? 

3\dz 

3x 

^n+l 

c 

0 

0 

1 

0 

Cn+1 

(E.  25) 
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and  the  fundamental  solution  matrix  may  be  written 


n(t)  = 


where  a,  b,  c  and  e  are  n  dimensional  vectors  and  d  is  a  scalar.  A  conjugate 
system  of  solutions  is  determined  by  its  initial  matrices,  now  taken  as 


'U(Q)' 

’i 

o' 

0 

1 

V<0) 

I 

0 

_0 

0_ 

(E.  27) 


for  when  the  initial  conditions  satisfy 
U'V  =  V'U 


(E.  28) 


then  so  do  the  solutions 


U  : 


:  U 
!  O 


i  r 

i  r 


TT  (t) 


!  V  ! 

-  O  '■ 


[VU{t)+A2{t)\  o 

[a(t)  +  c(t)j 

[n21(t)  +  Vt)l  pi  • 


1  ! 


(E. 29) 


J 


L- 


Then  the  determinant  of  U  is 

f  "l2«>]  0 


dei  U  =  det 


[a(t)  +  c<t)  ] 


3  det  [tTjjU)  +  n12(tj]  (E.  30) 
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rll(t) 

0  ' 

n12(t)  bit) 

( 

a'(t) 

1 

c'(t)  d{t> 

1T(  0)  =  I  . 

(E. 26) 

V2l{t) 

0  1 

i 

n22(t)  e(t) 

0 

0  1 

0  1 

i 

I 


JL 


The  reduced  system 


so  that  if  Hj^t)  +  tt^U) 


is  non-singular,  then  so  is  U. 


' 

- 

*i 

r 

UR 

n12(t) 

UR(0) 

UR(0) 

1  i 

«  ' 
> 

_ i 

n^it) 

n22^ 

J 

< 

o 

.vr(0>. 

"  L  ‘  _ 

(E. 3 1) 


is  a  conjugate  system  of  solutions  which  forms  a  field  for  the  reduced  problem 
i^[rll<t)  +  n12^t)J  is  non-singular.  Thus,  the  field  in  (n+l)-space  implies 
the  existence  of  the  field  in  (n+2)-space. 

The  system  (E.21)  results  when  V(0)  is  taken  as 

0 

0  0 


HE,  32) 


V(0) 


b  Z 


+  I 
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APPENDIX  F 

THE  LEAST  SINGULAR  SQUARE  MATRIX  OF  AN  nxq  MATRIX 


Let  n  be  greater  than  q,  and  let  each  row  vector  V.  of  the  matrix  be  normalized 
to  1, 


‘v'v 1  i 

i=  i 


V?.  -  1,  i  =  1,  ....  n. 


(F.  1) 


The  method  proceeds  by  construction,  so  let 


WK  =  *1V1  + 


+  aKVK 


(F.2) 


be  a  unit  vector  made  up  of  the  first  K  most  orthogonal  vectors,  and  find 
the  maximum  projection  (V  ,W„)  of  VV'  on  a  row  vector  V  .  m  >K.  The 
projection  is  the  cosine  of  the  angle  between  Vm  and  W^,  which,  when 
maximized,  gives  the  smallest  angle  between  Vm  and  the  K-dimensional 
Euclidian  space  determined  by  the  K  vectors  in  (F.  2),  Thus,  the  problem 
is  that  of  finding  the  minimum  of  the  maximum  projections  of  the  remaining 
(n-m)  row  vectors,  and  using  this  to  form  Wj^+j. 

Consider  the  maximization  problem  represented  by 


f  (a)  =  (V  ,WJ  =  maximum,  m  >  K, 
m  m  K 


(F .  3) 


subject  to  the  constraint 


(wK,wK)  =  1. 


Introduce  a  Lagrange  multiplier  X  and  form  the  expression 

P,a.M  »  <VWK)+*[1-(WK.WK)]. 


(F.  4) 


(F.5) 


1 
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or  more  simply 


b  =  $„  (2\  a), 
m  Is. 

(F .  9) 

from  which 

-1 

**  ■  [«k]  V 

(F.  10) 

Now  is  the  left-hand  side  of  (F.  8),  which,  when  multiplied  by  a,  gives  f  , 
so  in  view  of  (F.7), 


2  Mb  'a)  =  f  2  =  b  •  U  'M  bm. 
m  m  m  L  K  J  m 
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Equation  {F.  11)  shows  that  the  maximum  and  minimum  projections  have 
opposite  signs  and  allows  the  squared  value  of  the  projection  to  be  evaluated 
in  terms  of  known  quantities.  This  is  good  enough  for  the  computer  mechani¬ 
zation  of  the  sch  roe.  which  almost  suggests  itself. 


Consider  the  (n  x  n)  matrix 


(Vj.Vj)  .  .  .  (Vj.Vn) 


«w  •  •  •  <w 


(F.  12) 


which  is  symmetric  and  has  ones  down  the  main  diagonal.  At  the  Kth  stage, 
^his  may  be  partitioned  into  four  submatrices 


r*K  b' 

LB*  cj. 

so  consider  the  product 


B 

*_1  ° 

I  B 

B* 

C 

0  I 

»V‘  c_ 

(F.  13) 


(F.  14) 


The  (m-K)1*1  row  of  the  lower  left  submatrix  times  the  (m-K)**1  column  of  B 
gives  the  square  of  the  projection  (F.  11).  This  is  easily  computed  for  each  of 
the  remaining  vectors,  and  the  minimum  of  these  identifies  the  most  orthogonal 
vector  to  the  K-dimensional  subspace. 


Now  suppose  that  rows  and  columns  are  interchanged  in  the  last  matrix  of  (F.  14) 
to  bring  the  minimizing  row  and  column  next  to  the  identity  matrix.  This  retains 
the  symmetry  of  C.  Then  elementary  column  operations  may  be  used  to  obtain 
the  (K+l)  x  (K+l)  identity  matrix  for  the  next  step. 
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APPENDIX  G 


PHOOF  THAT  DISCONTINUITIES  IN  THE  MULTIPLIERS 
CAN  HE  DETERMINED  AT  EITHER  END  OK  A  CONSTRAINED  SUBARC 


The  geometry  is  shown  in  Figure  G-l.  It  is 

P 


to  he  shown  that  the  value 
p  ^  is  obtained  by  following 
either  path  A  or  path  B,  Path 
A  assumes  that  the  multipliers 
are  continuous  at  the  initial 
point  t,  and  discontinuous  at 
the  endpoint  t.,.  Path  B  assumes 
the  discontinuities  are  at  the 
initial  point. 


Figure  G  1 . 

Geometry  of  the  Problem 


The  necessary  conditions  for  the  constrained  subarc  are 


-P 

0 


f(x,  u) 


IM  *  Sf 

0  a  r 


I  V-  1 
t  3a  ; 


Jadl 

I-*  I  ( 
1  3X  1 


af 


,  af  ,  O  ,  3G  .  3H,  . 

P  au  +  air  +  P-aU  '  aCT (x*  u-  p-  **> 


(G.  1) 
(G.  2) 

(G.  3) 


0 


G  (x,  u) 


(G.  4) 


The  last  two  equations  are  "solved1  to  obtain 


u  =  u(x) 

,  -  p)  -  -  rjfe  <p'  I  +  >• 

1  9u  ) 

Substitution  of  (G.  5)  and  (G.  6)  into  (G.  1)  and  (G.  2)  gives  the  reduced 
differential  equations  for  the  constrained  subarc: 


x  =  f(x,  u(:<)) 


"9fo 

ax 

O 

r  iL.  /  ■  n 

3u 

(«-) 

+ 

_2L.  &L 

ax  aG  \  ax  1 

- 

51 r 

- 

L  sir  J 

Equation  (G.9)  has  been  added  for  convenience. 

Now  (G.  7)  and  (G.  8)  ar  e  "uncoupled"  since  p  does  not  appear  in 
(G.  7).  Since  x  =  x(t)  .and  tue  coefficients  of  p  in  (G.8)  contain 
x  alone,  (G.  8)  and  (G.  9)  form  a  linear  first-order  homogeneous 
system  of  differential  equations  with  time  varying  coefficients, 

ip  = 

where  \ 1/  is  the  n+i  dimensional  vector 


Let  <t  be  the  fundamental  solution  matrix  of  (G .  10)  with  $(tj)  =  I, 
tlie  (n+1)  x  (n  +  1 )  identity  matrix.  Then  the  value  of  computed 
along  path  A  is 


~vi 

=  <Mu) 

a 

-  V 

1  ax  1 

Ci 

i 

0 

(G.  12) 


Aloxig  path  B, 


—  — 

ion  /n 

Po 

i 

a 

"B 

_Po_ 

1  *<9 

L  1  J 

(G.  13) 


If  (G.  12)  and  (G.  13)  give  the  same  value  for  p2  it  follows  that 


l*>2 | 

i*m 

1  dx  | 

=  45  (t2> 

1  ax  j 

_ 

_ 

(G.  14) 


and  hence  that 


9G(t) 

9x 


is  a  solution  of  (G .  § )  v.'ith  p 


this,  consider  the  vector  identity 


0.  To  show 


3G 

bx 


(x,  u(x)) 


aG  +  aG  au 

ax  au  ax 


o 


(G.  15) 
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since 


9u  _  1  8G 

3  x  I  aG  i  ,'j  x 

I  au  I 


(G.  16) 


results  from  equations  (G.  1)  and  (G.  5).  On  the  other  hand. 


.  / 

1 9G1  - 

a  1 

[f '  ss\  - 

a2a  f  + 

af.'  ju' 

ie  r 

\  ax  |  “ 

ax  | 

l  Sxl 
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ax 

ax  ax 

iauJ. 

SC' 

ax  • 


(G.  17) 


Additionally, 


Combining  (G.  15)  -  (G.  18),  one  finds 


(G.  18, 


d_  9G' 
dt  ax 


af 
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af. 

aG 

o| 

ax 

'  (act 

Bu 

ax 

ax  ) 

m 

lau  1 

(G.  19) 


which  is  the  same  as  (G.  8)  with  p  =0  and  p  »  25.  Thus  — - 

*o  r  9x  ax 

is  a  solution  of  (G.  8)  with  pQ  =  0,  and  the  values  of  p0  computed 
along  paths  A  and  B  are  the  same. 
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